MODULE chem_sources

  USE grid_dims, ONLY: &
       maxgrds

  USE mem_chem1, ONLY: &
       nsrc,           &
       max_ntimes_src, &
       chem1_vars,     &       ! Type
       chem1_src_vars          ! Type

  USE mem_aer1, ONLY: &
       aer1_vars               ! Type

  USE mem_plume_chem1, ONLY: &
       plume_vars,           & ! Type
       plume_mean_vars         ! Type

  USE mod_chem_plumerise_scalar, ONLY: &
       plumerise_driver        ! Subroutine



  IMPLICIT NONE

  CHARACTER(len=256)                                  :: SRCMAPFN
  CHARACTER(len=32)                                   :: def_proc_src 
  ! def_proc_src='last_sources' !if does not exist => keep the current sources
  ! def_proc_src='stop' 	!if does not exist => stop the execution
  REAL                                                :: prfrq             
  INTEGER                                             :: plumerise 
  REAL                                                :: srctime1=0.
  REAL                                                :: srctime2=0.

  PRIVATE

  LOGICAL                                             :: got_init_cond = .FALSE.   
  LOGICAL                                             :: got_srcfiles_inv = .FALSE.
  INTEGER                                             :: src_swap=0
  INTEGER            , PARAMETER                      :: maxsrcfiles=1500
  CHARACTER(len=256) , DIMENSION(maxsrcfiles,maxgrds) :: fnames_src
  CHARACTER(len=14)  , DIMENSION(maxsrcfiles,maxgrds) :: itotdate_src
  REAL               , DIMENSION(maxsrcfiles,maxgrds) :: src_times
  INTEGER            , DIMENSION(maxgrds)             :: nsrcfiles
  INTEGER            , DIMENSION(maxgrds)             :: next_srcfile
  INTEGER            , DIMENSION(max_ntimes_src,nsrc) :: actual_time_index


  !- arrays for the diurnal cycle of emission
  TYPE cycle_emission
     DOUBLE PRECISION, POINTER, DIMENSION(:,:) :: dcnorma_inv,emission_rate
  END TYPE cycle_emission

  LOGICAL                                            :: emiss_cycle_alloc = .FALSE.
  TYPE(cycle_emission) , ALLOCATABLE, DIMENSION(:,:) :: emiss_cycle
  REAL                 , PARAMETER                   :: emiss_cycle_time = 86400.

  PUBLIC :: le_sourcemaps,   & ! Subroutine
            sources_driver,  & ! Subroutine
            sources,         & ! Subroutine
            SRCMAPFN,        &
            def_proc_src,    &
            prfrq,           &
            plumerise,       &
            srctime1,        &
            srctime2

CONTAINS

  !----------------------------------------------------------------------------

  SUBROUTINE le_sourcemaps(ng,m1,m2,m3,time,iyear1,imonth1,idate1,itime1,&
                           ngrids,timmax,chem_nspecies,spc_chem_alloc,   &
                           src,off,nsrc,nvert_src,chem1_src_g,bburn,     & 
                           spc_chem_name,on,chemical_mechanism,          &
                           emiss_ajust,co,aer_nspecies,spc_aer_alloc,    &
                           spc_aer_name,urban,nucle,accum,src_name,      &
                           chemistry,ntimes_src,aer1_g,nmodes,aerosol,   &
                           nveg_agreg,plume_mean_g,nzpmax,dzt,rtgt,      &
                           transport,aer_bburn,plume_g,tropical_forest,  &
                           boreal_forest,savannah,grassland,diur_cycle)

! (DMK) NOT USED
! USE mem_chem1
!
! (DMK) USED
! USE mem_grid, ONLY: &
!         time,          &
!         iyear1,        &
!         imonth1,       &
!         idate1,        &
!         itime1,        &
!         ngrids,        &
!         timmax
! (DMK) NOT USED
!    USE mem_grid, ONLY: &
!         grid_g,        &
!         runtype
!
! (DMK) NOT USED
!  USE mem_basic, ONLY: &
!         basic_g

    ! original
    INTEGER , INTENT(IN) :: ng
    INTEGER , INTENT(IN) :: m1
    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    
    ! grid_dims
    INTEGER, INTENT(IN) :: nzpmax

    ! mem_grid
    REAL    , INTENT(IN) :: time
    INTEGER , INTENT(IN) :: iyear1
    INTEGER , INTENT(IN) :: imonth1
    INTEGER , INTENT(IN) :: idate1
    INTEGER , INTENT(IN) :: itime1
    INTEGER , INTENT(IN) :: ngrids
    REAL    , INTENT(IN) :: timmax
    REAL    , INTENT(IN) :: dzt(nzpmax)
    REAL    , INTENT(IN) :: rtgt(m2,m3)

    ! chem1_list
    INTEGER          , INTENT(IN) :: chem_nspecies
    INTEGER          , INTENT(IN) :: spc_chem_alloc(6,chem_nspecies)
    INTEGER          , INTENT(IN) :: src
    CHARACTER(LEN=8) , INTENT(IN) :: spc_chem_name(chem_nspecies)
    INTEGER          , INTENT(IN) :: on
    INTEGER          , INTENT(IN) :: off
    CHARACTER(LEN=24), INTENT(IN) :: chemical_mechanism
    REAL             , INTENT(IN) :: emiss_ajust(chem_nspecies)
    INTEGER          , INTENT(IN) :: CO
    INTEGER          , INTENT(IN) :: transport

    ! mem_chem1
    INTEGER              , INTENT(IN)    :: nsrc
    INTEGER              , INTENT(IN)    :: nvert_src(nsrc,maxgrds)
    TYPE(chem1_src_vars) , INTENT(INOUT) :: chem1_src_g(2,nsrc,chem_nspecies,maxgrds)
    INTEGER              , INTENT(IN)    :: bburn
    CHARACTER(LEN=20)    , INTENT(IN)    :: src_name(nsrc)
    INTEGER              , INTENT(IN)    :: chemistry
    INTEGER              , INTENT(IN)    :: ntimes_src(nsrc)
    INTEGER              , INTENT(IN)    :: diur_cycle(nsrc)

    ! aer1_list
    INTEGER          , INTENT(IN) :: aer_nspecies
    INTEGER          , INTENT(IN) :: spc_aer_alloc(6,nmodes,aer_nspecies)
    CHARACTER(LEN=8) , INTENT(IN) :: spc_aer_name(nmodes,aer_nspecies)
    INTEGER          , INTENT(IN) :: urban
    INTEGER          , INTENT(IN) :: nucle
    INTEGER          , INTENT(IN) :: accum
    INTEGER          , INTENT(IN) :: aer_bburn

    ! mem_aer1
    TYPE (aer1_vars) , INTENT(INOUT) :: aer1_g(nmodes,aer_nspecies,maxgrds)
    INTEGER          , INTENT(IN)    :: nmodes
    INTEGER          , INTENT(IN)    :: aerosol

    ! mem_plume_chem1
    INTEGER               , INTENT(IN)    :: nveg_agreg
    TYPE (plume_mean_vars), INTENT(INOUT) :: plume_mean_g(nveg_agreg,maxgrds)
    TYPE (plume_vars)     , INTENT(INOUT) :: plume_g(nveg_agreg,chem_nspecies,maxgrds)
    INTEGER               , INTENT(IN)    :: tropical_forest
    INTEGER               , INTENT(IN)    :: boreal_forest
    INTEGER               , INTENT(IN)    :: savannah
    INTEGER               , INTENT(IN)    :: grassland


    !- local var
    INTEGER :: iunit,isrc,k1,k2,isrctime,iwtdo
    CHARACTER(len=256)  :: fname
    CHARACTER(len=2)   :: cgrid

    !    call makefnam(fname,SRCMAPFN//' ',time,iyear1,imonth1,idate1,itime1*100,'T',cgrid,'vfm')

    IF(SRCMAPFN(1:LEN_TRIM(SRCMAPFN)) == 'NONE' .OR. SRCMAPFN(1:LEN_TRIM(SRCMAPFN)) == 'none') RETURN

    IF( got_srcfiles_inv == .FALSE.) THEN    
       CALL src_file_inv (SRCMAPFN,iyear1,imonth1,idate1,itime1,ngrids,time,timmax,nsrc,diur_cycle)
!       CALL init_actual_time_index
       CALL init_actual_time_index(nsrc,ntimes_src)
    ENDIF

    ! def_proc_src='last_sources'  ! < move to namelist


    !- check if files required does exist and decide what to do in case not
    CALL check_src_files(ng,iwtdo)
    IF(iwtdo == 0) RETURN 


    !- loop at source files (isrctime will be > 1, only if will use linterp for any source)
    DO isrctime=1,MAXVAL(ntimes_src(1:nsrc))


       !- swap sources:  copy time level 2->1 and read the next data for the time level 2
       IF(src_swap == 1 .AND. MAXVAL(ntimes_src(1:nsrc)) > 1 .AND. isrctime==1) THEN 
!         CALL  swap_sources(ng,m1,m2,m3,time)
          CALL swap_sources(ng,m1,m2,m3,time,chem_nspecies,spc_chem_alloc, &
                            src,off,nsrc,nvert_src,chem1_src_g,bburn)

          next_srcfile(ng) = next_srcfile(ng) + 1 
          CYCLE
       ENDIF
       !- next file to open
       fname=fnames_src(next_srcfile(ng),ng)

       !- read emission dataset using V-format
       !if(trim(fname(len_trim(fname)-2:len_trim(fname))) == 'vfm') then

!       CALL read_sources_vfm(ng,m1,m2,m3,iyear1,imonth1,idate1,isrctime&
!            ,fname(1:LEN_TRIM(fname)))

       CALL read_sources_vfm(ng,m1,m2,m3,iyear1,imonth1,idate1,isrctime,    &
                             fname(1:LEN_TRIM(fname)),chem_nspecies,        &
                             spc_chem_alloc,spc_chem_name,src,on,off,       &
                             chemical_mechanism,emiss_ajust,co,aer_nspecies,&
                             spc_aer_alloc,spc_aer_name,urban,nucle,accum,  &
                             nsrc,chem1_src_g,src_name,chemistry,aer1_g,    &
                             nmodes,aerosol,nveg_agreg,plume_mean_g)

       !else;! read HDF;!endif


       !- split bburn emissions into flaming/smoldering parts
       IF(plumerise == 1) &
!            CALL emis_flam_smold(ng,m1,m2,m3,isrctime)
            CALL emis_flam_smold(ng,m1,m2,m3,isrctime,              &
                              nsrc,chem1_src_g,bburn,chem_nspecies, &
                              spc_chem_alloc,src,on,off,transport,  &
                              aer1_g,aerosol,aer_nspecies,          &
                              spc_aer_alloc,nmodes,aer_bburn,       &
                              plume_mean_g,plume_g,tropical_forest, &
                              boreal_forest,savannah,grassland,     &
                              nveg_agreg)


       !-----
       !- to use the mass conservation fix
       !- convert from [kg/m^2]to mixing ratio expressed in [ppbm = 1.e9 kg/kg]
       !call convert_to_mixing_ratio(ng,m1,m2,m3,isrctime)

       !- convert from [kg/m^2]to tracer density expressed in [ 1.e9 kg/m3]

!       CALL convert_to_tracer_density(ng,m1,m2,m3,isrctime)
       CALL convert_to_tracer_density(ng,m1,m2,m3,isrctime,             &
                                      nzpmax,dzt,rtgt,nsrc,             &
                                      chem1_src_g,bburn,chem_nspecies,  &
                                      spc_chem_alloc,src,off,transport, &
                                      aer1_g,aerosol,aer_nspecies,      &
                                      spc_aer_alloc,nmodes,aer_bburn) 
       !-----

       !- next time (don't change the position of this line)
       IF(isrctime==1) next_srcfile(ng) = next_srcfile(ng) + 1 

    ENDDO
    !- update srctime for the next time of reading/update    
    srctime1	 =   src_times(next_srcfile(ng)-1,1)
    srctime2	 =   src_times(next_srcfile(ng)  ,1)

    IF(ng==ngrids) src_swap = 1


    PRINT*,'next_srcfile,srctime1-2=',next_srcfile(ng),srctime1,srctime2
    CALL flush(6)


  END SUBROUTINE le_sourcemaps

  !----------------------------------------------------------------------
  SUBROUTINE src_file_inv(srcpref,iyear1,imonth1,idate1,itime1,ngrids, &
                          time,timmax,nsrc,diur_cycle)

    ! original
    CHARACTER(len=*) , INTENT(IN) :: srcpref
    INTEGER          , INTENT(IN) :: iyear1
    INTEGER          , INTENT(IN) :: imonth1
    INTEGER          , INTENT(IN) :: idate1
    INTEGER          , INTENT(IN) :: itime1
    INTEGER          , INTENT(IN) :: ngrids
    REAL             , INTENT(IN) :: time
    REAL             , INTENT(IN) :: timmax

    ! mem_chem1
    INTEGER , INTENT(IN) :: nsrc
    INTEGER , INTENT(IN) :: diur_cycle(nsrc)

    INTEGER :: nc,nf,lnf,nvftot, ng,it,isrc,nfstart
    INTEGER :: inyear,inmonth,indate,inhour


    CHARACTER(len=256), DIMENSION(maxsrcfiles) :: fnames
    CHARACTER(len=256) :: vpref
    CHARACTER(len=14)  :: itotdate,itotdate_current
    CHARACTER(len=2)   :: cgrid
    REAL(kind=8) :: secs_init,secs_src


    fnames(1:maxsrcfiles)= 'XXXXXXXXXXXXXXXX'
    ! Get abs seconds of run start

    CALL date_abs_secs2(iyear1,imonth1,idate1,itime1*100,secs_init)

    ! get the current time  
    CALL date_add_to(iyear1,imonth1,idate1,itime1*100  &
         ,time,'s',inyear,inmonth,indate,inhour)

    CALL date_make_big (inyear,inmonth,indate,inhour,itotdate_current)

    ! Go through src files and make inventory

    nc=LEN_TRIM(srcpref)
    nvftot=-1
    vpref=srcpref
    !print*,vpref; call flush(6)
    DO ng=1,ngrids

       WRITE(cgrid,'(a1,i1)') 'g',ng

       CALL RAMS_filelist(fnames,TRIM(vpref)  &
            //'*'//cgrid//'*.vfm',nvftot)

       IF(nvftot > maxsrcfiles) THEN
          STOP ' too many sources files'   
       ENDIF

       nsrcfiles(ng)=0
       DO nf=1,nvftot
          lnf=LEN_TRIM(fnames(nf))
          !print*,lnf,fnames(nf)

          READ(fnames(nf)(lnf-23:lnf-6),20) inyear,inmonth,indate,inhour
20        FORMAT(i4,1x,i2,1x,i2,1x,i6)

          ! print*,inyear,inmonth,indate,inhour
          CALL date_make_big(inyear,inmonth,indate,inhour,itotdate)

          nsrcfiles(ng)=nsrcfiles(ng)+1
          fnames_src(nsrcfiles(ng),ng)=fnames(nf)
          itotdate_src(nsrcfiles(ng),ng)=itotdate

          CALL date_abs_secs2(inyear,inmonth,indate,inhour,secs_src)
          src_times(nsrcfiles(ng),ng)=secs_src - secs_init

       ENDDO

       CALL RAMS_dintsort(nsrcfiles(ng),itotdate_src(:,ng),fnames_src(:,ng))

       !  start printing section
       !--------------------------------------------------------------

       PRINT*,' '
       PRINT*,' '
       PRINT*,' '
       PRINT*,'-------------------------------------------------------------'
       PRINT*,'-----------  Sources File Input Inventory --for --- GRID=', ng
       PRINT*,'-------------------------------------------------------------'
       DO nf=1,nsrcfiles(ng)
          PRINT 8,  nf, itotdate_src(nf,ng),src_times(nf,ng) ,TRIM(fnames_src(nf,ng))
       ENDDO
8      FORMAT(i4,1x,a16,1x,f10.0,2x,a)
       PRINT*,'------------------------------------------------------'

    ENDDO ! ngrids
    !call flush(6)

    IF(ngrids > 1 .AND. INT((SUM(nsrcfiles(1:ngrids)))/ nsrcfiles(1)) .NE. ngrids) THEN
       PRINT*,'The number of src files for each grid, MUST be the same:'
       PRINT*,'nsrcfiles(1:ngrids)=',nsrcfiles(1:ngrids)
       STOP 999
    ENDIF

    !- Are there enough src files available(for now only 1 grid is considered)
    !
    ng=1
    IF(src_times(nsrcfiles(ng),ng) < timmax) THEN 
       PRINT*,'=============================================================='
       PRINT*,'Warning: Not enough files for the entire time integration'
       PRINT*,'Warning: were found, model will continue.'
       PRINT*,'=============================================================='
    ENDIF

    !- perform some initializations (for now only 1 grid is considered)
    ng = 1 
    next_srcfile(1:ngrids) = 0
        !print*,'nvftot: ',nvftot
    DO nf=1,nvftot
        !print*,'todates: ',itotdate_src(nf,ng) , itotdate_current
       IF(itotdate_src(nf,ng) == itotdate_current) THEN
          next_srcfile(1:ngrids) = nf
          EXIT
       ELSEIF(itotdate_src(nf,ng) > itotdate_current) THEN
          next_srcfile(1:ngrids) = nf-1
          EXIT
       ENDIF
        !print*,'next_srcfile: ',next_srcfile(1:ngrids)
    ENDDO
        !print*,'next_srcfile: ',next_srcfile(1:ngrids)


!--(DMK-CCATT)---------------------------------------------------------
  IF (next_srcfile(ng) < 1) STOP ' next_srcfile < 1 '
  IF (next_srcfile(ng) > nvftot-1 .AND. SUM(diur_cycle(1:nsrc)) < 3 ) &
         STOP 'next_srcfile > nvftot-1'
!--(DMK-original)------------------------------------------------------
!    IF (next_srcfile(ng) < 1 .OR. next_srcfile(ng) > nvftot-1) &
!         STOP ' next_srcfile < 1 .or. next_srcfile > nvftot-1'
!--(DMK-CCATT-END)-----------------------------------------------------

    srctime1     =   src_times(next_srcfile(ng)  ,ng)
    srctime2     =   src_times(next_srcfile(ng)+1,ng)

!--(DMK-CCATT)---------------------------------------------------------
  IF (next_srcfile(ng) > nvftot-1 .AND. SUM(diur_cycle(1:nsrc)) == 3 ) &
   srctime2 = srctime1 + 86400. 
!--(DMK-CCATT-END)-----------------------------------------------------

    !- fill src_times arrays above nvftot with valid numbers (to use in 
    !- case of forecast or not more available data)
    DO ng=1,ngrids
       DO nf=nsrcfiles(ng)+1,maxsrcfiles
          src_times(nf,ng) = src_times(nf-1,ng) + (srctime2-srctime1)
       ENDDO
    ENDDO

    got_srcfiles_inv = .TRUE.

    ng=1
    print*,'next_srcfile=',next_srcfile(ng),srctime1,srctime2,itotdate_current
    print*,'next_srcfile= ',TRIM(fnames_src(next_srcfile(ng),ng))
    call flush(6)

  END SUBROUTINE src_file_inv
  !--------------------------------------------------------------

  SUBROUTINE init_actual_time_index(nsrc,ntimes_src)

    INTEGER , INTENT(IN) :: nsrc
    INTEGER , INTENT(IN) :: ntimes_src(nsrc)
    

    INTEGER :: it

    !- index to control memory access of src arrays, because some
    !- arrays have 2 time levels and others only 1 time level
    !- for time 1, the memory position is always allocated for all sources types
    it = 1
    actual_time_index(it,1:nsrc) = 1
    !- for the second, will depend on if linterp is wanted or not.
    !- if linterp is not desired for any sources, actual_time_index = 1
    !- and will have actually one only memory position 
    it = 2
    actual_time_index(it,1:nsrc) = ntimes_src(1:nsrc) 

  END SUBROUTINE init_actual_time_index
  !--------------------------------------------------------------

  SUBROUTINE check_src_files(ng,iwtdo)

    ! original
    INTEGER , INTENT(INOUT) :: iwtdo
    INTEGER , INTENT(IN)    :: ng

    iwtdo = 1

    !-check if is not greater the max number defined
    IF(next_srcfile(ng)+1 > maxsrcfiles) THEN
       PRINT*,'next_srcfile(ng)+1 > maxsrcfiles'
       CALL flush(6)
       STOP 
    ENDIF
    print*,' next_srcfile(ng) , nsrcfiles(ng)',next_srcfile(ng) , nsrcfiles(ng)
    call flush(6)

    IF(next_srcfile(ng) > nsrcfiles(ng)) THEN 
!    IF(next_srcfile(ng)+1 > nsrcfiles(ng)) THEN 

       !- situation 1 : stop model execution
       IF(def_proc_src(1:LEN_TRIM(def_proc_src)) == 'STOP' ) THEN
          PRINT*,'-----------------------------------------------------'
          PRINT*, 'Not src files available'
          PRINT*, 'stop at read_sources_vfm le_sourcemaps'
          CALL flush(6)
	  STOP

   !- situation 2 :keep the current sources
       ELSEIF(def_proc_src(1:LEN_TRIM(def_proc_src)) == 'LAST_SOURCES' )THEN

	  iwtdo=0

	  next_srcfile(ng)=next_srcfile(ng)+1
   !- update srctime for the next time of reading/update    
          srctime1	 =   src_times(next_srcfile(ng)-1,1)
          srctime2	 =   src_times(next_srcfile(ng)  ,1)

          !print*,'-----------------------------------------------------'
          PRINT*,'Not src files available:'
	  PRINT*,'using previous day sources, model will continue ...'
	  PRINT*,'srctime1 and 2=',srctime1,srctime2
	  CALL flush(6)


          RETURN
       ENDIF

    ENDIF


  END SUBROUTINE check_src_files
  !--------------------------------------------------------------

  SUBROUTINE swap_sources(ng,m1,m2,m3,time,chem_nspecies,spc_chem_alloc, &
                          src,off,nsrc,nvert_src,chem1_src_g,bburn)

! (DMK) NOT USED
!    USE chem1_list, ONLY:           &
!         spc_chem_name =>spc_name,  &
!         on,                        &
!         offline,                   &
!         transport
!
! (DMK) USED
!    USE chem1_list, ONLY:           &
!         chem_nspecies=>nspecies,   &
!         spc_chem_alloc=>spc_alloc, &
!         src,                       &
!         off
!
! (DMK) USED
!    USE mem_chem1, ONLY:               &
!         nvert_src=>chem1_src_z_dim_g, &
!         chem1_src_g,                  &
!         bburn
!
! (DMK) NOT USED
!    USE mem_chem1, ONLY:               &
!         bioge,                        &
!         antro                         &
!         chem1_g,                      &
!         nsrc,                         &
!
! (DMK) NOT USED
!    USE mem_aer1, ONLY:                   &
!         aer_nvert_src=>aer1_src_z_dim_g, &
!         aer1_g,                          &
!         AEROSOL
!
! (DMK) NOT USED
!    USE aer1_list, ONLY:           &
!         aer_nspecies=>nspecies,   &
!         spc_aer_alloc=>spc_alloc, &
!         nmodes,                   &
!         spc_aer_name =>aer_name,  &
!         aer_bburn => bburn,       &
!         aer_sdust => sdust,       &
!         aer_urban => urban,       &
!         aer_bioge => bioge,       &
!         aer_marin => marin 

    ! original
    INTEGER , INTENT(IN) :: ng
    INTEGER , INTENT(IN) :: m1
    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    REAL    , INTENT(IN) :: time

    ! chem1_list
    INTEGER , INTENT(IN) :: chem_nspecies
    INTEGER , INTENT(IN) :: spc_chem_alloc(6,chem_nspecies)
    INTEGER , INTENT(IN) :: src
    INTEGER , INTENT(IN) :: off

    ! mem_chem1
    INTEGER             , INTENT(IN)    :: nsrc
    INTEGER             , INTENT(IN)    :: nvert_src(nsrc,maxgrds)
    TYPE(chem1_src_vars), INTENT(INOUT) :: chem1_src_g(2,nsrc,chem_nspecies,maxgrds)
    INTEGER             , INTENT(IN)    :: bburn

    INTEGER ispc,isrc
    INTEGER, PARAMETER :: it1=1 & ! time level 1
         ,it2=2   ! time level 2

    DO isrc=1,nsrc  

       !- if the time level 2 uses the same memory allocation area
       !- of time level 1 => nothing to do
       IF(actual_time_index(it2,isrc) == 1) CYCLE
       IF(isrc == bburn) STOP 444

       !- else: make the swap     
       DO ispc=1,chem_nspecies
          IF(spc_chem_alloc(src,ispc) == off) CYCLE

          chem1_src_g(it1,isrc,ispc,ng)%sc_src(1:nvert_src(isrc,ng),1:m2,1:m3) = &
               chem1_src_g(it2,isrc,ispc,ng)%sc_src(1:nvert_src(isrc,ng),1:m2,1:m3) 
       ENDDO
    ENDDO

    !- aerosol section still need to be done

    PRINT*,'--> source swapped done at time (h)=',time/3600.


  END SUBROUTINE swap_sources

  !--------------------------------------------------------------

  SUBROUTINE read_sources_vfm(ng,m1,m2,m3,iyear,imon,iday,isrctime,fname, &
                              chem_nspecies,spc_chem_alloc,spc_chem_name, &
                              src,on,off,chemical_mechanism,emiss_ajust,  &
                              co,aer_nspecies,spc_aer_alloc,spc_aer_name, &
                              urban,nucle,accum,nsrc,chem1_src_g,src_name,&
                              chemistry,aer1_g,nmodes,aerosol,nveg_agreg, &
                              plume_mean_g) 

! (DMK) NOT USED
!    USE chem1_list, ONLY:           &
!         offline,                   &
!
! (DMK) USED
!    USE chem1_list, ONLY:           &
!         chem_nspecies=>nspecies,   &
!         spc_chem_alloc=>spc_alloc, &
!         spc_chem_name =>spc_name,  &
!         src,                       &
!         on,                        &
!         off,                       &
!         chemical_mechanism,        &
!         emiss_ajust,               &
!         CO
!
! (DMK) NOT USED
!    USE aer1_list, ONLY:           &
!         nmodes,                   &
!         sdust,                    &
!         bburn,                    &
!         bioge,                    &
!         marin,                    &
!         coarse,                   &
!
! (DMK) USED
!    USE aer1_list, ONLY:           &
!         aer_nspecies=>nspecies,   &
!         spc_aer_alloc=>spc_alloc, & 
!         spc_aer_name =>aer_name,  &
!         urban,                    &
!         nucle,                    &
!         accum
!
! (DMK) USED
!    USE mem_chem1, ONLY: &
!         chem1_src_g,    &
!         src_name,       &
!         chemistry
!
! (DMK) USED
!    USE mem_aer1, ONLY:  &
!         aer1_g,         &
!         nmodes,         &
!         aerosol
!
! (DMK) USED
!    USE mem_plume_chem1, ONLY: &
!         plume_mean_g

    ! original
    INTEGER       , INTENT(IN) :: ng
    INTEGER       , INTENT(IN) :: m1
    INTEGER       , INTENT(IN) :: m2
    INTEGER       , INTENT(IN) :: m3
    INTEGER       , INTENT(IN) :: iyear
    INTEGER       , INTENT(IN) :: imon
    INTEGER       , INTENT(IN) :: iday
    INTEGER       , INTENT(IN) :: isrctime
    CHARACTER*(*) , INTENT(IN) :: fname

    ! chem1_list
    INTEGER          , INTENT(IN) :: chem_nspecies
    INTEGER          , INTENT(IN) :: spc_chem_alloc(6,chem_nspecies)
    CHARACTER(LEN=8) , INTENT(IN) :: spc_chem_name(chem_nspecies)
    INTEGER          , INTENT(IN) :: src
    INTEGER          , INTENT(IN) :: on
    INTEGER          , INTENT(IN) :: off
    CHARACTER(LEN=24), INTENT(IN) :: chemical_mechanism
    REAL             , INTENT(IN) :: emiss_ajust(chem_nspecies)
    INTEGER          , INTENT(IN) :: CO

    ! aer1_list
    INTEGER          , INTENT(IN) :: aer_nspecies
    INTEGER          , INTENT(IN) :: spc_aer_alloc(6,nmodes,aer_nspecies)
    CHARACTER(LEN=8) , INTENT(IN) :: spc_aer_name(nmodes,aer_nspecies)
    INTEGER          , INTENT(IN) :: urban
    INTEGER          , INTENT(IN) :: nucle
    INTEGER          , INTENT(IN) :: accum

    ! mem_chem1
    INTEGER              , INTENT(IN)    :: nsrc
    TYPE(chem1_src_vars) , INTENT(INOUT) :: chem1_src_g(2,nsrc,chem_nspecies,maxgrds)
    CHARACTER(LEN=20)    , INTENT(IN)    :: src_name(nsrc)
    INTEGER              , INTENT(IN)    :: chemistry

    ! mem_aer1
    TYPE (aer1_vars) , INTENT(INOUT) :: aer1_g(nmodes,aer_nspecies,maxgrds)
    INTEGER          , INTENT(IN)    :: nmodes
    INTEGER          , INTENT(IN)    :: aerosol

    ! mem_plume_chem1
    INTEGER               , INTENT(IN)    :: nveg_agreg
    TYPE (plume_mean_vars), INTENT(INOUT) :: plume_mean_g(nveg_agreg,maxgrds)

    !- local var
    INTEGER :: iunit,ispc,isrc,iveg_agreg,nvert,ihour,itim
    CHARACTER(len=20) read_spc_name,read_src_name,section,read_units,date
    CHARACTER(len=20) read_mean,read_veg_name
    INTEGER read_ident_chem_mec,read_ident_src,dummy,read_ident_veg
    INTEGER read_ident_aer,read_aer_mode, imode
    INTEGER nspecies,nxp,nyp,i,j
    REAL dep_glon(2), dep_glat(2)
    REAL, ALLOCATABLE, DIMENSION(:,:) :: src_dummy_2d
    CHARACTER(len=32) :: chemical_mechanism_test
    LOGICAL :: there

    PRINT*,'-----------------------------------------------------'
    PRINT*,'Opening source file= ',fname(1:LEN_TRIM(fname))

    INQUIRE(file=fname,exist=there)
    IF(.NOT.there) THEN
       PRINT*,'-----------------------------------------------------'
       PRINT*,'FILE= ',TRIM(fname(1:LEN_TRIM(fname))),' does not exist.'

       PRINT*, 'STOP at read_sources_vfm le_sourcemaps'
       PRINT*,'-----------------------------------------------------'
       CALL flush(6)
       STOP	   
    ENDIF

    !- initial attributions/allocations
    nspecies = chem_nspecies + aer_nspecies*nmodes

    ihour = 0
    nvert = 1
    iunit = 2
    ALLOCATE (src_dummy_2d(m2,m3));src_dummy_2d=0.


    !-need to zerout aerosol sources, before reading the emission of the new day
    IF(AEROSOL == on) THEN
       DO ispc=1,aer_nspecies
          DO imode=1,nmodes
             IF(spc_aer_alloc(src,imode,ispc) == on) aer1_g(imode,ispc,ng)%sc_src = 0.
          ENDDO
       ENDDO
    ENDIF

    !- open the source file
    OPEN(UNIT=iunit,FILE=fname,FORM='formatted',STATUS='old') 
    ispc = 0 ; isrc = 0; iveg_agreg = 0

    READ(iunit,*)  nxp,(dep_glon(i),i=1,2)
    !print*,nxp,(dep_glon(i),i=1,2)
    !call flush(6)
    READ(iunit,*) nyp,(dep_glat(i),i=1,2)
    !print*, nyp,(dep_glat(i),i=1,2)
    !call flush(6)

    READ(iunit,*)  date 
    !- test if the source data is for the chemical mechanism that will be used:
    READ(iunit,*)  chemical_mechanism_test
    IF(TRIM( chemical_mechanism_test ) /=  TRIM(chemical_mechanism)) THEN
       PRINT*,'chemical mechanism not the same as expected'
       PRINT*,'expected=',TRIM(chemical_mechanism(1:LEN_TRIM(chemical_mechanism)))
       PRINT*,'read    =',TRIM(chemical_mechanism_test(1:LEN_TRIM(chemical_mechanism_test)))
       STOP 'wrong chem mechanism at chem_sources'
    ELSE
       PRINT*,'Reading sources for chem mechanism= ',TRIM(chemical_mechanism(1:LEN_TRIM(chemical_mechanism)))
    ENDIF

    DO i=1,nspecies*nsrc*5

       READ(iunit,*,END=100) section

       !- emission section  --------------------------------
       IF(TRIM(section) == 'chemistry') THEN
          READ(iunit,*)   read_spc_name &
               , read_ident_chem_mec &
               , read_src_name       &
               , read_ident_src      &
               , read_units

          ispc = read_ident_chem_mec
          isrc = read_ident_src
          itim = actual_time_index(isrctime,isrc)

!!!!!if(spc_chem_alloc(src,ispc) == off) then
          IF(.NOT. ASSOCIATED( chem1_src_g(itim,isrc,ispc,ng)%sc_src)) THEN
             PRINT*, 'CHEM memory not allocated for specie=',ispc,TRIM(read_spc_name)
             PRINT*, 'at time=',actual_time_index(isrctime,isrc)
             PRINT*, 'stop at read_sources_vfm routine'
             STOP 100
          ENDIF

          PRINT*,'ispc=',ispc,TRIM(spc_chem_name(ispc))
          PRINT*,'reading emission of ', TRIM(spc_chem_name(ispc)), ' source: ', TRIM(src_name(isrc))

          CALL vfirec(iunit,chem1_src_g(itim,isrc,ispc,ng)%sc_src(1,:,:),m2*m3,'LIN')


          !- this is not correct for reading fluxes (for example, CO2 assimilation)
          !where (chem1_src_g(itim,isrc,ispc,ng)%sc_src(1,:,:) < 0.) &
          !       chem1_src_g(itim,isrc,ispc,ng)%sc_src(1,:,:) = 0.


          !- Biogenic CO only for tracer runs
!          IF(CHEMISTRY > 0 .AND.  ispc == CO .AND. spc_chem_alloc(src,CO) == ON ) THEN           
          IF(CHEMISTRY > 0 .AND.  ispc == CO .AND. spc_chem_alloc(3,CO) == ON ) THEN           
             !print*,'bioge=',bioge; call flush(6)
             !chem1_src_g(itim,bioge,CO,ng)%sc_src(1,:,:) = 0.
             !srf- cheque erro compilador intel
             chem1_src_g(itim,3,CO,ng)%sc_src(1,:,:) = 0.
          ENDIF
          !
          !- ajust emissions 
          chem1_src_g(itim,isrc,ispc,ng)%sc_src(1,:,:) = emiss_ajust(ispc)*&
               chem1_src_g(itim,isrc,ispc,ng)%sc_src(1,:,:)


!ajustando emissoes antropicas
!        if(TRIM(src_name(isrc))=='antro') chem1_src_g(itim,isrc,ispc,ng)%sc_src(1,:,:)= &
!	                           0.3333*chem1_src_g(itim,isrc,ispc,ng)%sc_src(1,:,:)
!-------------
!
          !- aerosol section  --------------------------------
       ELSEIF(TRIM(section) == 'aerosol' .AND. AEROSOL == on) THEN
          READ(iunit,*)   read_spc_name       &
               , read_ident_aer      &
               , read_aer_mode       &
               , read_src_name       &
               , read_ident_src      &
               , read_units
          ispc = read_ident_aer
          imode= read_aer_mode
          isrc = read_ident_src
          PRINT*,'ispc=',ispc,imode,TRIM(spc_aer_name(imode,ispc))

          IF(spc_aer_alloc(src,imode,ispc) == off)STOP ' AER memory not allocated for '

          PRINT*,'reading emission of ', TRIM(spc_aer_name(imode,ispc)), ' source: ', TRIM(src_name(isrc))

          CALL vfirec(iunit,src_dummy_2d,m2*m3,'LIN')

          !- this is not correct for reading fluxes (for example, CO2 assimilation)
          !where (src_dummy_2d(:,:) < 0.) src_dummy_2d(:,:) =0.

          aer1_g(imode,ispc,ng)%sc_src(1,:,:) = aer1_g(imode,ispc,ng)%sc_src(1,:,:) +src_dummy_2d(:,:)



          !- plume section  --------------------------------
       ELSEIF(TRIM(section) == 'plume' .AND. plumerise == on) THEN

          READ(iunit,*)   read_mean &
               , dummy &
               , read_veg_name       &
               , read_ident_veg      &
               , read_units
          PRINT*,'reading ',read_mean,' for ',read_veg_name

          iveg_agreg = read_ident_veg       
          IF(TRIM(read_mean) == 'mean_fct' ) THEN
             CALL vfirec(iunit,plume_mean_g(iveg_agreg,ng)%flam_frac(:,:),m2*m3,'LIN')
             WHERE(plume_mean_g(iveg_agreg,ng)%flam_frac(:,:) < 0.) &
                  plume_mean_g(iveg_agreg,ng)%flam_frac(:,:) = 0.

          ELSEIF(TRIM(read_mean) == 'firesize' ) THEN
             CALL vfirec(iunit,plume_mean_g(iveg_agreg,ng)%fire_size(:,:),m2*m3,'LIN')
             WHERE(plume_mean_g(iveg_agreg,ng)%fire_size(:,:) < 0.) &
                  plume_mean_g(iveg_agreg,ng)%fire_size(:,:) = 0.

          ELSE
             STOP ' model not yet prepared for flam-frac of each one of species'
          ENDIF
       ELSEIF(TRIM(section) /= 'plume' .AND. plumerise == on)  THEN

          STOP 'unknow error at routine read_sources_vfm'
       ENDIF

    ENDDO
100 CONTINUE
    PRINT*,'-----------------------------------------------------'
    CLOSE (iunit) 

    IF(AEROSOL == on) THEN
       !- special section for aerosol sulfate 
       !- mode Aitken = 50% mode Accumulation
       !* ref: Stier et al., The aerosol-climate model ECHAM5-HAM. Atmos.
       !  Chem. Phys., 5,1125-1156,2005.
       IF(spc_aer_alloc(src,nucle,urban) == on) &
            aer1_g(nucle,urban,ng)%sc_src(1,:,:)=0.5*aer1_g(accum,urban,ng)%sc_src(1,:,:)

       IF(spc_aer_alloc(src,accum,urban) == on) &
            aer1_g(accum,urban,ng)%sc_src(1,:,:)=0.5*aer1_g(accum,urban,ng)%sc_src(1,:,:)

       !print*,'aer1=',aer1_g(accum,urban,ng)%sc_src(1,:,:)

    ENDIF

    ! do j=1,m3
    !  do i=1,m2
    !    print*,'frac=',plume_mean_g(:,ng)%flam_frac(i,j)
    !    print*,'size=',plume_mean_g(:,ng)%fire_size(i,j)
    ! enddo;enddo
    ! stop 222
    DEALLOCATE (src_dummy_2d)
  END SUBROUTINE read_sources_vfm

  !-------------------------------------------------------------
  ! (DMK) NOT USED
  !
  ! subroutine convert_to_mixing_ratio(ng,m1,m2,m3,isrctime) 
  !  use mem_basic, only: basic_g
  !  use mem_grid, only: dzt,grid_g
  !  use mem_chem1 
  !  use chem1_list, only : chem_nspecies=>nspecies &
  !                        ,spc_chem_alloc=>spc_alloc&
  !		          ,spc_chem_name =>spc_name,src,on,off,offline&
  !			  ,transport
  !
  !  use mem_aer1
  !  use aer1_list, only : aer_nspecies=>nspecies &
  !                       ,spc_aer_alloc=>spc_alloc, nmodes&
  !                       ,spc_aer_name =>aer_name &
  !			 ,aer_bburn => bburn 
  !  implicit none
  !
  !  integer,intent(IN) :: ng,m1,m2,m3,isrctime
  !  
  !  
  !  ! local var
  !  !Fator de conversao de unidades    
  !  real,parameter :: fcu =1.e+9 !=> ppbm 
  !  !LPCE
  !  real :: rhodz_inv,dz
  !  integer :: i,j,ksrc,isrc,ispc,imode,itim
  !  
  !  ksrc=2 !surface level of emission in the model
  !  
  !  do j=1,m3
  !     do i=1,m2
  !        
  !        ! Todas as unidades estao em kg/m2/dia => use 'dz' em vez de 'vol'
  !        !    vol = 1./(dxt(i,j)*dyt(i,j)*dzt(k))*rtgt(i,j)
  !
  !        dz        = grid_g(ng)%rtgt(i,j)/dzt(ksrc) ! dzt=1/(z(k)-z(k-1))	
  !	  rhodz_inv = 1./(dz*basic_g(ng)%dn0(ksrc,i,j))
  !	  do ispc=1,chem_nspecies
  !           if(spc_chem_alloc(src,ispc) == off) cycle
  !
  !                !- convert from kg/m^2  to  density (kg[gas]/m^3*1.e9)
  ! 	     	  do isrc=1,nsrc
  !
  !		      itim = actual_time_index(isrctime,isrc)
  !
  !	              chem1_src_g(itim,isrc,ispc,ng)%sc_src(1,i,j) = &
  !		      chem1_src_g(itim,isrc,ispc,ng)%sc_src(1,i,j) * fcu * rhodz_inv
  !                enddo
  !                ! copy smoldering emission from  bburn to ksrc
  !		  itim = actual_time_index(isrctime,bburn)
  !                chem1_src_g(itim,bburn,ispc,ng)%sc_src(ksrc,i,j) = &
  !		  chem1_src_g(itim,bburn,ispc,ng)%sc_src(1   ,i,j)
  !
  !	  enddo
  !	enddo 
  !  enddo
  !  
  ! !- aerosol section 
  ! IF(AEROSOL > 0 ) then
  !  do j=1,m3
  !     do i=1,m2
  !        
  !        ! Todas as unidades estao em kg/m2/dia => use 'dz' em vez de 'vol'
  !        !    vol = 1./(dxt(i,j)*dyt(i,j)*dzt(k))*rtgt(i,j)
  !
  !        dz        = grid_g(ng)%rtgt(i,j)/dzt(ksrc) ! dzt=1/(z(k)-z(k-1))	
  !	  rhodz_inv = 1./(dz*basic_g(ng)%dn0(ksrc,i,j))
  ! 	  do ispc=1,aer_nspecies
  !	    do imode=1,nmodes
  !	     
  !	      if(spc_aer_alloc(src,imode,ispc)       == off .or. &
  !               spc_aer_alloc(transport,imode,ispc) == off) cycle
  !            
  !	      !- convert from kg/m^2  to  density (kg[aer]/m^3*1.e9)
  !           
  !	     
  !	      aer1_g(imode,ispc,ng)%sc_src(1,i,j) = &
  !	      aer1_g(imode,ispc,ng)%sc_src(1,i,j) * fcu * rhodz_inv
  !          
  !          enddo
  !        enddo
  !
  !	 
  !	  ! copy smoldering emission from aerosol bburn to ksrc
  !	  do imode=1,nmodes
  !	     
  !	     if(spc_aer_alloc(src,imode,aer_bburn)       == off .or. &
  !               spc_aer_alloc(transport,imode,aer_bburn) == off) cycle
  !         
  !	     aer1_g(imode,aer_bburn,ng)%sc_src(ksrc,i,j) = &
  !	     aer1_g(imode,aer_bburn,ng)%sc_src(1   ,i,j)
  !        enddo
  !
  !	enddo 
  !  enddo
  ! ENDIF
  !
  ! end subroutine convert_to_mixing_ratio
  !-------------------------------------------------------------

  SUBROUTINE convert_to_tracer_density(ng,m1,m2,m3,isrctime,               &
                                       nzpmax,dzt,rtgt,nsrc,chem1_src_g,   &
                                       bburn,chem_nspecies,spc_chem_alloc, &
                                       src,off,transport,aer1_g,aerosol,   &
                                       aer_nspecies,spc_aer_alloc,nmodes,  &
                                       aer_bburn) 

! (DMK) NOT USED
!
!    USE mem_basic, ONLY: &
!         basic_g
!
! (DMK) USED
!    USE mem_grid, ONLY: &
!         dzt,           &
!         grid_g
!
! DMK) USED
!    USE mem_chem1, ONLY: &
!         chem1_src_g,    &
!         bburn
!
! (DMK) NOT USED
!    USE chem1_list, ONLY:           &
!         spc_chem_name =>spc_name,  &
!         on,                        &
!         offline,                   &
!
! (DMK) USED
!    USE chem1_list, ONLY:           &
!         chem_nspecies=>nspecies,   &
!         spc_chem_alloc=>spc_alloc, &
!         src,                       &
!         off,                       &
!         transport
!
! (DMK) USED
!    USE mem_aer1, ONLY: &
!         aer1_g,        &
!         aerosol
!
! (DMK) NOT USED
!    USE aer1_list, ONLY:           &
!         spc_aer_name =>aer_name,  &
!
! (DMK) USED
!    USE aer1_list, ONLY:           &
!         aer_nspecies=>nspecies,   &
!         spc_aer_alloc=>spc_alloc, &
!         nmodes,                   &
!         aer_bburn => bburn 

    ! original
    INTEGER , INTENT(IN) :: ng
    INTEGER , INTENT(IN) :: m1
    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    INTEGER , INTENT(IN) :: isrctime

    ! grid_dims
    INTEGER , INTENT(IN) :: nzpmax

    ! mem_grid
    REAL    , INTENT(IN) :: dzt(nzpmax)
    REAL    , INTENT(IN) :: rtgt(m2,m3)

    ! chem1_list
    INTEGER , INTENT(IN) :: chem_nspecies
    INTEGER , INTENT(IN) :: spc_chem_alloc(6,chem_nspecies)
    INTEGER , INTENT(IN) :: src
    INTEGER , INTENT(IN) :: off
    INTEGER , INTENT(IN) :: transport

    ! aer1_list
    INTEGER , INTENT(IN) :: aer_nspecies
    INTEGER , INTENT(IN) :: spc_aer_alloc(6,nmodes,aer_nspecies)
    INTEGER , INTENT(IN) :: nmodes
    INTEGER , INTENT(IN) :: aer_bburn

    ! mem_chem1
    INTEGER              , INTENT(IN)    :: nsrc
    TYPE(chem1_src_vars) , INTENT(INOUT) :: chem1_src_g(2,nsrc,chem_nspecies,maxgrds)
    INTEGER              , INTENT(IN)    :: bburn

    ! mem_aer1
    TYPE (aer1_vars) , INTENT(INOUT) :: aer1_g(nmodes,aer_nspecies,maxgrds)
    INTEGER          , INTENT(IN)    :: aerosol

    ! local var
    !Fator de conversao de unidades    
    REAL,PARAMETER :: fcu =1.e+9 !=> ppbm 
    !LPCE
    REAL :: dz_inv,dz
    INTEGER :: i,j,ksrc,isrc,ispc,imode,itim

    ksrc=2 !surface level of emission in the model

    !- chemistry section

    DO j=1,m3
       DO i=1,m2

          ! Todas as unidades estao em kg/m2/dia => use 'dz' em vez de 'vol'
          !    vol = 1./(dxt(i,j)*dyt(i,j)*dzt(k))*rtgt(i,j)

!          dz        = grid_g(ng)%rtgt(i,j)/dzt(ksrc) ! dzt=1/(z(k)-z(k-1))	
          dz        = rtgt(i,j)/dzt(ksrc) ! dzt=1/(z(k)-z(k-1))	
	  dz_inv    = 1./dz
	  DO ispc=1,chem_nspecies
             IF(spc_chem_alloc(src,ispc) == off) CYCLE

             !- convert from kg/m^2  to  density (kg[gas]/m^3*1.e9)
             DO isrc=1,nsrc
                itim = actual_time_index(isrctime,isrc)

                chem1_src_g(itim,isrc,ispc,ng)%sc_src(1,i,j) = &
                     chem1_src_g(itim,isrc,ispc,ng)%sc_src(1,i,j) * fcu * dz_inv
             ENDDO
             ! copy smoldering emission from  bburn to ksrc
             itim = actual_time_index(isrctime,bburn)

             chem1_src_g(itim,bburn,ispc,ng)%sc_src(ksrc,i,j) = &
		  chem1_src_g(itim,bburn,ispc,ng)%sc_src(1   ,i,j)

	  ENDDO
       ENDDO
    ENDDO

    !- aerosol section 
    IF(AEROSOL > 0 ) THEN
       DO j=1,m3
          DO i=1,m2

             ! Todas as unidades estao em kg/m2/dia => use 'dz' em vez de 'vol'
             !    vol = 1./(dxt(i,j)*dyt(i,j)*dzt(k))*rtgt(i,j)

!             dz        = grid_g(ng)%rtgt(i,j)/dzt(ksrc) ! dzt=1/(z(k)-z(k-1))	
             dz        = rtgt(i,j)/dzt(ksrc) ! dzt=1/(z(k)-z(k-1))	
             dz_inv    = 1./dz
             DO ispc=1,aer_nspecies
                DO imode=1,nmodes

                   IF(spc_aer_alloc(src,imode,ispc)       == off .OR. &
                        spc_aer_alloc(transport,imode,ispc) == off) CYCLE

                   !- convert from kg/m^2  to  density (kg[aer]/m^3*1.e9)


                   aer1_g(imode,ispc,ng)%sc_src(1,i,j) = &
                        aer1_g(imode,ispc,ng)%sc_src(1,i,j) * fcu * dz_inv

                ENDDO
             ENDDO


             ! copy smoldering emission from aerosol bburn to ksrc
             DO imode=1,nmodes

                IF(spc_aer_alloc(src,imode,aer_bburn)       == off .OR. &
                     spc_aer_alloc(transport,imode,aer_bburn) == off) CYCLE

                aer1_g(imode,aer_bburn,ng)%sc_src(ksrc,i,j) = &
                     aer1_g(imode,aer_bburn,ng)%sc_src(1   ,i,j)
             ENDDO

          ENDDO
       ENDDO
    ENDIF

  END SUBROUTINE convert_to_tracer_density

  !----------------------------------------------------------------------

  SUBROUTINE  emis_flam_smold(ng,n1,n2,n3,isrctime,                 &
                              nsrc,chem1_src_g,bburn,chem_nspecies, &
                              spc_chem_alloc,src,on,off,transport,  &
                              aer1_g,aerosol,aer_nspecies,          &
                              spc_aer_alloc,nmodes,aer_bburn,       &
                              plume_mean_g,plume_g,tropical_forest, &
                              boreal_forest,savannah,grassland,     &
                              nveg_agreg)

! (DMK) USED    
!    USE mem_chem1, ONLY: &
!         chem1_src_g,    &
!         bburn
!
! (DMK) NOT USED
!    USE chem1_list, ONLY:           &
!         spc_chem_name =>spc_name,  &
!         offline,                   &
!
! (DMK) USED
!    USE chem1_list, ONLY:           &
!         chem_nspecies=>nspecies,   &
!         spc_chem_alloc=>spc_alloc, &
!         src,                       &
!         on,                        &
!         off,                       &
!         transport
!
! (DMK) USED
!    USE mem_aer1, ONLY: &
!         aer1_g,        &
!         aerosol
!
! (DMK) NOT USED
!    USE aer1_list, ONLY:           &
!         aer_nspecies=>nspecies,   &
!         spc_aer_name =>aer_name,  &

! (DMK) USED
!    USE aer1_list, ONLY:           &
!         spc_aer_alloc=>spc_alloc, &
!         nmodes,                   &
!         aer_bburn => bburn ! aerosol type
!
! (DMK) USED
!    USE mem_plume_chem1, ONLY: &
!         plume_mean_g,         &
!         plume_g,              &
!         tropical_forest,      &
!         boreal_forest,        &
!         savannah,             &
!         grassland,            &
!         nveg_agreg

    ! original
    INTEGER , INTENT(IN) :: ng
    INTEGER , INTENT(IN) :: n1
    INTEGER , INTENT(IN) :: n2
    INTEGER , INTENT(IN) :: n3
    INTEGER , INTENT(IN) :: isrctime

    ! chem1_list
    INTEGER          , INTENT(IN) :: chem_nspecies
    INTEGER          , INTENT(IN) :: spc_chem_alloc(6,chem_nspecies)
    INTEGER          , INTENT(IN) :: src
    INTEGER          , INTENT(IN) :: on
    INTEGER          , INTENT(IN) :: off
    INTEGER          , INTENT(IN) :: transport

    ! aer1_list
    INTEGER          , INTENT(IN) :: nmodes
    INTEGER          , INTENT(IN) :: aer_nspecies
    INTEGER          , INTENT(IN) :: spc_aer_alloc(6,nmodes,aer_nspecies)
    INTEGER          , INTENT(IN) :: aer_bburn

    ! mem_chem1
    INTEGER              , INTENT(IN)    :: nsrc
    TYPE(chem1_src_vars) , INTENT(INOUT) :: chem1_src_g(2,nsrc,chem_nspecies,maxgrds)
    INTEGER              , INTENT(IN)    :: bburn


    ! mem_aer1
    TYPE (aer1_vars) , INTENT(INOUT) :: aer1_g(nmodes,aer_nspecies,maxgrds)
    INTEGER          , INTENT(IN)    :: aerosol

    ! mem_plume_chem1
    INTEGER               , INTENT(IN)    :: nveg_agreg
    TYPE (plume_mean_vars), INTENT(INOUT) :: plume_mean_g(nveg_agreg,maxgrds)
    TYPE (plume_vars)     , INTENT(INOUT) :: plume_g(nveg_agreg,chem_nspecies,maxgrds)
    INTEGER               , INTENT(IN)    :: tropical_forest
    INTEGER               , INTENT(IN)    :: boreal_forest
    INTEGER               , INTENT(IN)    :: savannah
    INTEGER               , INTENT(IN)    :: grassland

    REAL,DIMENSION(n2,n3) :: smold_frac 
    INTEGER iv,ispc,i,j,imode,itim
    INTEGER:: imean_plume
    imean_plume = 1 !change this at alloc_plume_chem1 routine also


    !- time index of memory allocation position        
    itim = actual_time_index(isrctime,bburn)
    IF(itim > 1) STOP'time level 2 not allowed when plumerise is used -emis_flam_smold'


    IF(imean_plume == on) THEN
       !-----  
       !- calcula a emissao smoldering e fatores para obtencao da fracao
       !- flaming em funcao da emissao smoldering

       smold_frac(:,:) = 1.- ( plume_mean_g(tropical_forest,ng)%flam_frac(:,:) + &
            plume_mean_g(boreal_forest  ,ng)%flam_frac(:,:) + &
            plume_mean_g(savannah	  ,ng)%flam_frac(:,:) + &
            plume_mean_g(grassland	  ,ng)%flam_frac(:,:)	)    


       !- chemistry section (only for bburn source)
       DO ispc = 1,chem_nspecies

  	  IF(spc_chem_alloc(src,ispc) /= on) CYCLE 

     !- convert from 'total' emisson to 'smoldering'
  	  chem1_src_g(itim,bburn,ispc,ng)%sc_src(1,:,:) = smold_frac(:,:) * &  
               chem1_src_g(itim,bburn,ispc,ng)%sc_src(1,:,:)
       ENDDO
       IF(AEROSOL > 0 ) THEN
          !- aerosol section (only for bburn aerosols)
          DO imode=1,nmodes

             IF(spc_aer_alloc(src,imode,aer_bburn)       == off .OR. &
                  spc_aer_alloc(transport,imode,aer_bburn) == off) CYCLE

             !- convert from 'total' emisson to 'smoldering'
             aer1_g(imode,aer_bburn,ng)%sc_src(1,:,:) = smold_frac(:,:) * & 
                  aer1_g(imode,aer_bburn,ng)%sc_src(1,:,:)

          ENDDO
       ENDIF

       !- convert from flaming fraction to relationship with phase smoldering emission
       DO iv = 1, nveg_agreg
          plume_mean_g(iv,ng)%flam_frac(:,:) = plume_mean_g(iv,ng)%flam_frac(:,:)/ &
               (1.e-8+smold_frac(:,:))
          !- flamming emission =  plume_mean_g(iv,ng)%flam_frac(:,:) * &
          !					 chem1_src_g(itim,bburn,ispc,ng)%sc_src(1,:,:)
       ENDDO

       !-----

    ELSE

       !-----      case where each specie has his own flaming fraction ----------------
       STOP 'aerosol emission not ready for this option'
       DO ispc = 1,chem_nspecies
          IF(spc_chem_alloc(src,ispc) /= on) CYCLE 
          smold_frac(:,:) = 1.- ( plume_g(tropical_forest,ispc,ng)%fct(:,:) + &
               plume_g(boreal_forest  ,ispc,ng)%fct(:,:) + &
               plume_g(savannah	  ,ispc,ng)%fct(:,:) + &
               plume_g(grassland	  ,ispc,ng)%fct(:,:)   )



          !- convert from 'total' emisson to 'smoldering'
          chem1_src_g(itim,bburn,ispc,ng)%sc_src(1,:,:) = smold_frac(:,:) * &  
               
               chem1_src_g(itim,bburn,ispc,ng)%sc_src(1,:,:)

          !- convert from flaming fraction to relationship with phase smoldering emission
          DO iv = 1, nveg_agreg
             plume_g(iv,ispc,ng)%fct(:,:) = plume_g(iv,ispc,ng)%fct(:,:)/ &
                  (1.e-8+smold_frac(:,:))
             !- flamming emission =  plume_g(iv,iscp,ng)%fct(:,:) * &
             !				      chem1_src_g(itim,bburn,ispc,ng)%sc_src(1,:,:)
          ENDDO

       ENDDO
       !-----
    ENDIF
  END SUBROUTINE  emis_flam_smold
  !----------------------------------------------------------------------

  SUBROUTINE sources_driver(mzp,mxp,myp,ia,iz,ja,jz, &
                            ngrid,mmxp,mmyp,maxgrds,cosz,theta,pp,pi0,      &
                            rv,dn0,up,vp,                                   &!srf-AWE
			    time,iyear1,imonth1,idate1,itime1,       &
                            dtlt,rtgt,lpw,glat,glon,zt,zm,dzt,ngrids,       &
                            nzpmax,pi180,&!nkp,ntime,zt_plumegen,zm_plumegen, &
                            !dzt_plumegen,dzm,dz,ncall,nrec,n_setgrid,       &
                            chem_nspecies,spc_chem_alloc,spc_chem_name,     &
                            src,on,off,transport,aer_nspecies,spc_aer_alloc,&
                            nmodes,aer_bburn,aer_sdust,aer_urban,aer_bioge, &
                            aer_marin,nvert_src,chem1_g,nsrc,bburn,         &
                            chem1_src_g,antro,bioge,chemistry,ntimes_src,   &
                            diur_cycle,aerosol,aer1_g,aer_nvert_src,        &
                            nveg_agreg,tropical_forest,boreal_forest,       &
                            savannah,grassland,plume_mean_g,dnp,iexev)

    ! original
    INTEGER , INTENT(IN) :: mzp
    INTEGER , INTENT(IN) :: mxp
    INTEGER , INTENT(IN) :: myp
    INTEGER , INTENT(IN) :: ia
    INTEGER , INTENT(IN) :: iz
    INTEGER , INTENT(IN) :: ja
    INTEGER , INTENT(IN) :: jz

    ! mem_all
    INTEGER , INTENT(IN) :: ngrid

    ! node_mod
    INTEGER , INTENT(IN) :: mmxp(maxgrds)
    INTEGER , INTENT(IN) :: mmyp(maxgrds)
    INTEGER , INTENT(IN) :: maxgrds

    ! mem_radiate
    REAL , INTENT(IN) :: cosz(mxp,myp)

    ! mem_basic
    REAL,    INTENT(IN) :: theta(mzp,mxp,myp)
    REAL,    INTENT(IN) :: pp(mzp,mxp,myp)
    REAL,    INTENT(IN) :: pi0(mzp,mxp,myp)
    REAL,    INTENT(IN) :: rv(mzp,mxp,myp)
    REAL,    INTENT(IN) :: up(mzp,mxp,myp)!srf-AWE
    REAL,    INTENT(IN) :: vp(mzp,mxp,myp)!srf-AWE
    REAL,    POINTER    :: dn0(:,:,:)

    ! mem_grid
    REAL    , INTENT(IN) :: time
    INTEGER , INTENT(IN) :: iyear1
    INTEGER , INTENT(IN) :: imonth1
    INTEGER , INTENT(IN) :: idate1
    INTEGER , INTENT(IN) :: itime1
    REAL    , INTENT(IN) :: dtlt
    REAL    , INTENT(IN) :: rtgt(mxp,myp)
    INTEGER , INTENT(IN) :: lpw(mxp,myp)
    REAL    , INTENT(IN) :: glat(mxp,myp)
    REAL    , INTENT(IN) :: glon(mxp,myp)
    REAL    , INTENT(IN) :: zt(nzpmax)
    REAL    , INTENT(IN) :: zm(nzpmax)
    REAL    , INTENT(IN) :: dzt(nzpmax)
    INTEGER , INTENT(IN) :: ngrids

    ! grid_dims
    INTEGER, INTENT(IN) :: nzpmax

    ! rconstants
    REAL , INTENT(IN) :: pi180

    ! plumegen_coms
    !INTEGER , INTENT(IN)    :: nkp 
    !INTEGER , INTENT(IN)    :: ntime 
    !REAL    , INTENT(INOUT) :: zt_plumegen(nkp) 
    !REAL    , INTENT(INOUT) :: zm_plumegen(nkp) 
    !REAL    , INTENT(INOUT) :: dzt_plumegen(nkp)
    !REAL    , INTENT(INOUT) :: dzm(nkp)
    !REAL    , INTENT(INOUT) :: dz 

    ! save attributes
    !INTEGER , INTENT(INOUT) :: ncall
    !INTEGER , INTENT(INOUT) :: nrec
    !INTEGER , INTENT(INOUT) :: n_setgrid

    ! chem1_list
    INTEGER          , INTENT(IN) :: chem_nspecies
    INTEGER          , INTENT(IN) :: spc_chem_alloc(6,chem_nspecies)
    INTEGER          , INTENT(IN) :: src
    CHARACTER(LEN=8) , INTENT(IN) :: spc_chem_name(chem_nspecies)
    INTEGER          , INTENT(IN) :: on
    INTEGER          , INTENT(IN) :: off
    INTEGER          , INTENT(IN) :: transport

    ! aer1_list
    INTEGER , INTENT(IN) :: aer_nspecies
    INTEGER , INTENT(IN) :: nmodes
    INTEGER , INTENT(IN) :: spc_aer_alloc(6,nmodes,aer_nspecies)
    INTEGER , INTENT(IN) :: aer_bburn
    INTEGER , INTENT(IN) :: aer_sdust
    INTEGER , INTENT(IN) :: aer_urban
    INTEGER , INTENT(IN) :: aer_bioge
    INTEGER , INTENT(IN) :: aer_marin

    ! mem_chem1
    INTEGER              , INTENT(IN)    :: nvert_src(nsrc,maxgrds)
    TYPE(chem1_vars)     , INTENT(INOUT) :: chem1_g(chem_nspecies,maxgrds)
    INTEGER              , INTENT(IN)    :: nsrc
    INTEGER              , INTENT(IN)    :: bburn
    TYPE(chem1_src_vars) , INTENT(INOUT) :: chem1_src_g(2,nsrc,chem_nspecies,maxgrds)
    INTEGER              , INTENT(IN)    :: antro
    INTEGER              , INTENT(IN)    :: bioge
    INTEGER              , INTENT(IN)    :: chemistry
    INTEGER              , INTENT(IN)    :: ntimes_src(nsrc)
    INTEGER              , INTENT(IN)    :: diur_cycle(nsrc)

    ! mem_aer1
    INTEGER          , INTENT(IN)    :: aerosol
    TYPE (aer1_vars) , INTENT(INOUT) :: aer1_g(nmodes,aer_nspecies,maxgrds)
    INTEGER          , INTENT(IN)    :: aer_nvert_src(aer_nspecies,maxgrds)

    ! mem_plume_chem1
    INTEGER               , INTENT(IN)    :: nveg_agreg
    INTEGER               , INTENT(IN)    :: tropical_forest
    INTEGER               , INTENT(IN)    :: boreal_forest
    INTEGER               , INTENT(IN)    :: savannah
    INTEGER               , INTENT(IN)    :: grassland
    TYPE (plume_mean_vars), INTENT(INOUT) :: plume_mean_g(nveg_agreg,maxgrds)

    ! mem_stilt
    REAL    , POINTER    :: dnp(:,:,:) ! in
    INTEGER , INTENT(IN) :: iexev

    IF( CHEMISTRY < 0 .OR. SRCMAPFN(1:LEN_TRIM(SRCMAPFN)) == 'NONE' &
         .OR. SRCMAPFN(1:LEN_TRIM(SRCMAPFN)) == 'none') RETURN

    !- allocation for diurnal cycle of emission arrays   and
    !- actual_time_index on nodes
    IF(.NOT. emiss_cycle_alloc) THEN
!       CALL alloc_emiss_cycle
       CALL alloc_emiss_cycle(mmxp,mmyp,ngrids,nsrc,maxgrds)
!       CALL init_actual_time_index
       CALL init_actual_time_index(nsrc,ntimes_src)
    ENDIF

    IF(MOD(time+0.001,emiss_cycle_time) .LT. dtlt .OR. time .LT. .01) & 
         CALL get_diurnal_cycle_normalized(ngrid,mxp,myp,ia,iz,ja,jz,dtlt,glat, &
                                           glon,imonth1,idate1,iyear1,itime1,   &
                                           antro,bioge,pi180)


    !-  call plumerise and sources  
    IF( plumerise == 1 ) &
         CALL plumerise_driver(mzp,mxp,myp,ia,iz,ja,jz,nzpmax,theta,pp,pi0,rv,up,vp, &!srf-AWE
	                       rtgt,lpw,     &
                               dtlt,time,ngrid,zt,zm,dzt,prfrq,srctime1,chem_nspecies,      &
                               spc_chem_alloc,src,on,off,nmodes,aer_nspecies,spc_aer_alloc, &
                               aer_bburn,maxgrds,nsrc,chem1_src_g,bburn,aer1_g,AEROSOL,     &
                               plume_mean_g,nveg_agreg,tropical_forest,boreal_forest,       &
                               savannah,grassland)!,nkp,ntime,zt_plumegen,dz,zm_plumegen,dzm, &
                               !dzt_plumegen,ncall,nrec,n_setgrid)

!    CALL sources(ngrid,mzp,mxp,myp,ia,iz,ja,jz,itime1,time,imonth1,idate1,iyear1 &
!         ,grid_g(ngrid)%glon)

    CALL sources(ngrid,mzp,mxp,myp,ia,iz,ja,jz,itime1,time,imonth1,idate1,iyear1,     &
                 glon,chem_nspecies,spc_chem_alloc,src,on,off,transport,nsrc,         &
                 nvert_src,chem1_src_g,chem1_g,bburn,bioge,antro,diur_cycle,aer_nvert_src,aer1_g,&
                 aerosol,aer_nspecies,spc_aer_alloc,nmodes,aer_bburn,aer_sdust,       &
                 aer_urban,aer_bioge,aer_marin,dnp,iexev,dn0,cosz)


  END SUBROUTINE  sources_driver

  !----------------------------------------------------------------------

  SUBROUTINE sources(ng,m1,m2,m3,ia,iz,ja,jz,itime1,time,imonth1,idate1,iyear1,glon,    &
                     chem_nspecies,spc_chem_alloc,src,on,off,transport,nsrc,nvert_src,  &
                     chem1_src_g,chem1_g,bburn,bioge,antro,diur_cycle,aer_nvert_src,    &
                     aer1_g,aerosol,aer_nspecies,spc_aer_alloc,nmodes,aer_bburn,        &
                     aer_sdust,aer_urban,aer_bioge,aer_marin,dnp,iexev,dn0,cosz)

! (DMK) NOT USED
!    USE chem1_list, ONLY :          &
!         spc_chem_name =>spc_name,  &
!         offline,                   &
!
! (DMK) USED
!    USE chem1_list, ONLY :          &
!         chem_nspecies=>nspecies,   &
!         spc_chem_alloc=>spc_alloc, &
!         src,                       &
!         on,                        &
!         off,                       &
!         transport
!
! (DMK) NOT USED
!    USE mem_chem1, ONLY: &
!         nsrc,                         &
!
! (DMK) USED
!    USE mem_chem1, ONLY: &
!         nvert_src=>chem1_src_z_dim_g, &
!         chem1_src_g,                  &
!         chem1_g,                      &
!         bburn,                        &
!         bioge,                        &
!         antro
!
! (DMK) USED
!    USE mem_aer1, ONLY:                   &
!         aer_nvert_src=>aer1_src_z_dim_g, &
!         aer1_g,                          &
!         AEROSOL
!
! (DMK) NOT USED
!    USE aer1_list, ONLY :          &
!        spc_aer_name =>aer_name,  &
!
! (DMK) USED
!    USE aer1_list, ONLY :          &
!         aer_nspecies=>nspecies,   &
!         spc_aer_alloc=>spc_alloc, &
!         nmodes,                   &
!         aer_bburn => bburn,       &
!         aer_sdust => sdust,       &
!         aer_urban => urban,       &
!         aer_bioge => bioge,       &
!         aer_marin => marin 
!
! (DMK) USED
!    USE mem_stilt, ONLY:  &
!         stilt_g,         &
!         iexev
!
! (DMK) USED
!    USE mem_basic, ONLY: &
!         basic_g
!
! (DMK) USED
!    USE mem_radiate, ONLY: &
!         radiate_g
!
! (DMK) NOT USED
!    USE mem_grid, ONLY: &
!         dtlt,          &
!         grid_g
!
! (DMK) NOT USED
!    USE node_mod, ONLY: &
!         mynum

    ! original
    INTEGER , INTENT(IN) :: ng
    INTEGER , INTENT(IN) :: m1
    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    INTEGER , INTENT(IN) :: ia
    INTEGER , INTENT(IN) :: iz
    INTEGER , INTENT(IN) :: ja
    INTEGER , INTENT(IN) :: jz
    REAL    , INTENT(IN) :: time
    INTEGER , INTENT(IN) :: imonth1
    INTEGER , INTENT(IN) :: idate1
    INTEGER , INTENT(IN) :: iyear1
    INTEGER , INTENT(IN) :: itime1

    ! chem1_list
    INTEGER , INTENT(IN) :: chem_nspecies
    INTEGER , INTENT(IN) :: spc_chem_alloc(6,chem_nspecies)
    INTEGER , INTENT(IN) :: src
    INTEGER , INTENT(IN) :: on
    INTEGER , INTENT(IN) :: off
    INTEGER , INTENT(IN) :: transport

    ! aer1_list
    INTEGER , INTENT(IN) :: aer_nspecies
    INTEGER , INTENT(IN) :: spc_aer_alloc(6,nmodes,aer_nspecies)
    INTEGER , INTENT(IN) :: nmodes
    INTEGER , INTENT(IN) :: aer_bburn
    INTEGER , INTENT(IN) :: aer_sdust
    INTEGER , INTENT(IN) :: aer_urban
    INTEGER , INTENT(IN) :: aer_bioge
    INTEGER , INTENT(IN) :: aer_marin

    ! mem_chem1
    INTEGER              , INTENT(IN)    :: nsrc
    INTEGER              , INTENT(IN)    :: nvert_src(nsrc,maxgrds)
    TYPE(chem1_src_vars) , INTENT(INOUT) :: chem1_src_g(2,nsrc,chem_nspecies,maxgrds)
    TYPE(chem1_vars)     , INTENT(INOUT) :: chem1_g(chem_nspecies,maxgrds)
    INTEGER              , INTENT(IN)    :: bburn
    INTEGER              , INTENT(IN)    :: bioge
    INTEGER              , INTENT(IN)    :: antro
    INTEGER              , INTENT(IN)    :: diur_cycle(nsrc)

    ! mem_aer1
    INTEGER         , INTENT(IN)    :: aer_nvert_src(aer_nspecies,maxgrds)
    TYPE(aer1_vars) , INTENT(INOUT) :: aer1_g(nmodes,aer_nspecies,maxgrds)
    INTEGER         , INTENT(IN)    :: aerosol

    ! mem_stilt
    REAL    , POINTER    :: dnp(:,:,:) ! in
    INTEGER , INTENT(IN) :: iexev

    ! mem_basic
    REAL , POINTER :: dn0(:,:,:) ! in

    ! mem_radiate
    REAL , INTENT(IN) :: cosz(m2,m3)

    INTEGER :: k_src,k_tend,it1,it2

    INTEGER,EXTERNAL :: julday

    DOUBLE PRECISION :: tlinterp
    !MFA
    REAL :: timeq2,timeq3,gglon,fuso,alfa(nsrc)
    INTEGER :: iweek,idays,j,i,k,ispc,isrc,k2,imode
    REAL :: tign,strtim,timeq,r_q,r_antro,real_time,jd
    REAL, DIMENSION(7) :: week_CYCLE
    !                     dia da semana:    SEG   TER	QUA   QUI   SEX   SAB  DOM  
    !                            iweek=     1      2	  3	4     5     6	 7
    !- dados cetesb/campinas/2005
    DATA (week_CYCLE(iweek),iweek=1,7) /1.1, 1.1, 1.1, 1.1, 1.1, 0.83,0.67/ !total = 7

    REAL rt(nsrc),rt_aer(aer_nspecies)
    REAL, PARAMETER :: bx_bburn  = 18.041288 * 3600., & !- pico em 18 UTC
         cx        = 2.5 * 3600.,       & ! 2.184936 * 3600., &
         rinti     = 0.8941851* 2.1813936e-8    , & ! 1/integral
         ax        = 2000.6038        , &
         bx_antro  = 9. *3600.  ,&  ! local time of peak 1 
         cx_antro  = 16.*3600.  ,&  ! local time of peak 2
         rsum      = 1.3848466E-05  !2.1311416E-05   ! 1/integral

    REAL, POINTER, SAVE,DIMENSION(:,:,:) :: rho_air
    REAL, DIMENSION(m2,m3) :: glon

    REAL local_cosz(m2,m3)

    !- nocturnal/background/constant emission for biogenic/urban-industrial-transp processes
    REAL, PARAMETER :: f_nct=0.15                  &! 15% per day
         , f_nct_dvd86400=f_nct/86400. &
         , um86400=1./86400.

    !- if using mass conservation fix : air dens changes with the time evolution
!    IF( iexev == 2 )  rho_air => stilt_g(ng)%dnp(:,:,:) 
    IF( iexev == 2 )  rho_air => dnp(:,:,:) 

    !- if not, air dens = air dens of basic state and  need to define for each when 
    !- have nested grids
!    IF( iexev == 1  ) rho_air => basic_g(ng)%dn0(:,:,:) 
    IF( iexev == 1  ) rho_air => dn0(:,:,:) 
    !
    !number of days of simulation
    idays = INT(( float(itime1)/100. + time/3600.)/24.+.00001)  

    IF(diur_cycle(bburn) == ON) THEN
       !-------------biomass burning diurnal cycle --------------------
       tign  = REAL(idays)*24.*3600.

       ! Modulacao da queimada media durante o ciclo diurno(unidade: 1/s)
       ! com a int( r_q dt) (0 - 24h)= 1.
       timeq= ( time + float(itime1)*0.01*3600. - tign )

       r_q  = rinti*( ax * EXP( -(timeq-bx_bburn)**2/(2.*cx**2) ) + 100. -  &
            5.6712963e-4*( timeq ))

       emiss_cycle(bburn,ng)%emission_rate(:,:)=r_q
       rt(bburn)= r_q
       alfa(bburn) = 0.
    ELSE
       emiss_cycle(bburn,ng)%emission_rate(:,:)= um86400 ! = 1./86400.
       rt(bburn)= um86400 ! = 1./86400.
       alfa(bburn) = 1.
    ENDIF


    IF(diur_cycle(antro) == ON) THEN

       !------------- anthropogenic diurnal cycle (industrial,urban, ...)
       real_time = time + float(itime1)*0.01*3600. !UTC

       ! weekly cycle
       ! week day
       !v1
       iweek= INT(((float(julday(imonth1,idate1,iyear1))/7. - &
            INT(julday(imonth1,idate1,iyear1)/7))*7.)) + 3
       IF(iweek.GT.7) iweek = iweek-7
       !v2
       !call Greg2Jul(0, idate1, imonth1, iyear1, jd)
       !jd=jd+float(nint(time/86400.))
       !iweek = int(AMOD(jd, 7.)) 
       !if(iweek < 1 .or. iweek > 7) stop 315

       !- diurnal cycle
       DO j=ja,jz
    	  DO i=ia,iz

             gglon = glon(i,j)
             fuso = INT(gglon/15)
             !-to better keep continuity
             idays = INT(( real_time +fuso*3600. )/86400.+.00001)  
             tign  = REAL(idays)*86400.

             timeq2= (real_time  -tign + fuso*3600.) - bx_antro
             timeq3= (real_time  -tign + fuso*3600.) - cx_antro
             !r_antro  =(exp(-((timeq2)**2)/((3.*3600) **2))+(exp(-((timeq3)**2)/((3.*3600) **2))+0.1))*rsum
             r_antro=(EXP(-(timeq2**2)/(18400.**2))+EXP(-(timeq3**2)/(18500.**2))+0.1)*rsum

             !- weekly + diurnal cycle
             r_antro = r_antro * week_CYCLE(iweek)

             emiss_cycle(antro,ng)%emission_rate(i,j)=r_antro 
             alfa(antro) = 0.
    	  ENDDO
       ENDDO
    ELSE
       !----------- sources linearly time interpolated
       IF(srctime2<=srctime1) STOP 'srctime2<=srctime1! linterp'
       tlinterp=DBLE(time-srctime1)/DBLE(srctime2-srctime1)
       emiss_cycle(antro,ng)%emission_rate(:,:)=tlinterp
       alfa(antro)= 1.
    ENDIF

    IF(diur_cycle(bioge) == ON) THEN

       !---------- sources with diurnal cycle and spatial dependence
       ! - using zenital angle from radiate routine/including constant background emission (f_cnt%)
!       emiss_cycle(bioge,ng)%emission_rate(:,:)= f_nct_dvd86400 + emiss_cycle(bioge,ng)%dcnorma_inv(:,:)&    
!            *MAX(0.,radiate_g(ng)%cosz(:,:)) * (1.-f_nct) 
       emiss_cycle(bioge,ng)%emission_rate(:,:)= f_nct_dvd86400 + emiss_cycle(bioge,ng)%dcnorma_inv(:,:)&    
            *MAX(0.,cosz(:,:)) * (1.-f_nct) 
       alfa(bioge) = 0.
    ELSE
       !----------- sources linearly time interpolated    
       IF(srctime2<=srctime1) STOP 'srctime2==srctime1! linterp'
       tlinterp=DBLE(time-srctime1)/DBLE(srctime2-srctime1)
       emiss_cycle(bioge,ng)%emission_rate(:,:) = tlinterp
       alfa(bioge) = 1.
    ENDIF


    ! print*,'emiss_cycle(antro,ng)%emission_rate=',emiss_cycle(antro,ng)%emission_rate(ia,ja)&
    ! ,srctime1,srctime2,tlinterp

    !-------------------------- perform emissions
    !- chemistry section 
    DO ispc=1,chem_nspecies

       IF(spc_chem_alloc(src,ispc) == off .OR. spc_chem_alloc(transport,ispc) == off)  CYCLE

       DO isrc=1,nsrc

          !- memory position of source array for each time level
          it1 = 1 ! always 1
          it2 = actual_time_index(2,isrc)!might be 1 or 2, will be 1 if linterp = OFF
          !=>chem1_src_g(it1,:,:,:)% = chem1_src_g(it2,:,:,:)%

          !print*,' nodes nvert=',isrc,nvert_src(isrc,ng)


          k_src = 1;	   k2 = 2     ! control for 2-dim src emission field
          IF(isrc == bburn )  k2 = m1-1  ! control for 3-dim src emission field

          DO k=2,k2

             IF(isrc == bburn ) k_src=k     ! control for 3-dim source emission field

             CALL source_to_tend_cycle ( m1,m2,m3,ia,iz,ja,jz   &
                  ,chem1_src_g(it1,isrc,ispc,ng)%sc_src &! source data at time level 1
                  ,chem1_src_g(it2,isrc,ispc,ng)%sc_src &! source data at time level 2
                  ,chem1_g    (	     ispc,ng)%sc_t   &! tendency array
                  ,nvert_src(isrc,ng)		     &! vertical size of source array
                  ,emiss_cycle(isrc,ng)%emission_rate   &! diurnal cycle of emission
                  ,alfa(isrc)                           &! alfa cte 
                  ,k_src				     &! vertical level of source array (where data is stored)
                  ,k				     &! vertical level of tendency
                  ,rho_air)			      ! air density (to convert to mixing ratio tendency)

          ENDDO
       ENDDO
    ENDDO

    !- aerosol section
    IF(AEROSOL > 0)  THEN
       !- still need implementation of the emission cycle with space dependence
       rt_aer(aer_bburn) = rt(bburn)
       rt_aer(aer_sdust) = 1.157407e-5
       rt_aer(aer_urban) = 1.157407e-5 !rt(antro)  ! < must be fixed later with actual diurnal cycle
       rt_aer(aer_bioge) = 1.157407e-5 
       rt_aer(aer_marin) = 1.157407e-5 

       DO ispc=1,aer_nspecies

          DO imode=1,nmodes

             IF(spc_aer_alloc(src      ,imode,ispc) == off .OR. &
                  spc_aer_alloc(transport,imode,ispc) == off) CYCLE


             k_src = 1;		  k2 = 2     ! control for 2-dim aerosol (bioge, antro, marin, sdust)
             IF(ispc == aer_bburn )  k2 = m1-1  ! control for 3-dim bburn aerosol

             DO k=2,k2

                IF(ispc == aer_bburn ) k_src=k     ! control  for 3-dim bburn aerosol

                CALL source_to_tend ( m1,m2,m3,ia,iz,ja,jz   &
                     ,aer1_g(imode,ispc,ng)%sc_src         	  &! source data
                     ,aer1_g(imode,ispc,ng)%sc_t		  &! tendency array
                     ,aer_nvert_src(ispc,ng)		  &! vertical size of source array
                     ,rt_aer(ispc)				  &! diurnal cycle of emission
                     ,k_src 				  &! vertical level of source array (where data is stored)
                     ,k					  &! vertical level of tendency
                     ,rho_air)				   ! air density (to convert to mixing ratio tendency)

             ENDDO
          ENDDO
       ENDDO
    ENDIF

  END SUBROUTINE sources
  !------------------------------------------------------------------

  SUBROUTINE source_to_tend(m1,m2,m3,ia,iz,ja,jz,sc_src,sc_t,nvert,rt,k_src,k_tend,rho_air)
   
    ! original
    INTEGER , INTENT(IN)    :: m1
    INTEGER,  INTENT(IN)    :: m2
    INTEGER,  INTENT(IN)    :: m3
    INTEGER,  INTENT(IN)    :: ia
    INTEGER,  INTENT(IN)    :: iz
    INTEGER,  INTENT(IN)    :: ja
    INTEGER,  INTENT(IN)    :: jz
    REAL,     INTENT(IN)    :: sc_src(nvert,m2,m3)
    REAL,     INTENT(INOUT) :: sc_t(m1,m2,m3)
    INTEGER,  INTENT(IN)    :: nvert
    REAL,     INTENT(IN)    :: rt
    INTEGER,  INTENT(IN)    :: k_src
    INTEGER,  INTENT(IN)    :: k_tend
    REAL,     INTENT(IN)    :: rho_air(m1,m2,m3)

    !- also this routine assumes that the source term is expressed in terms of  density (kg/m3)
    !- and, then, the conversion to mixing ratio (division by air density) is needed.
    !- k_src express the level where the src data is stored in the sc_src array
    sc_t(k_tend,ia:iz,ja:jz) = sc_t(k_tend,ia:iz,ja:jz) +  sc_src(k_src,ia:iz,ja:jz)*rt /  &
                                !- the level of air density
                                !- must be the same of tendency
                               rho_air(k_tend,ia:iz,ja:jz) 

  END SUBROUTINE source_to_tend


  !------------------------------------------------------------------

  SUBROUTINE source_to_tend_cycle(m1,m2,m3,ia,iz,ja,jz,sc_src1,sc_src2,sc_t,nvert, &
                                  rt,alfa,k_src,k_tend,rho_air)

    ! original
    INTEGER          , INTENT(IN)    :: m1
    INTEGER          , INTENT(IN)    :: m2
    INTEGER          , INTENT(IN)    :: m3
    INTEGER          , INTENT(IN)    :: ia
    INTEGER          , INTENT(IN)    :: iz
    INTEGER          , INTENT(IN)    :: ja
    INTEGER          , INTENT(IN)    :: jz
    REAL             , INTENT(IN)    :: sc_src1(nvert,m2,m3)
    REAL             , INTENT(IN)    :: sc_src2(nvert,m2,m3)
    REAL             , INTENT(INOUT) :: sc_t(m1,m2,m3)
    INTEGER          , INTENT(IN)    :: nvert
    DOUBLE PRECISION , INTENT(IN)    :: rt(m2,m3)
    REAL             , INTENT(IN)    :: alfa
    INTEGER          , INTENT(IN)    :: k_src
    INTEGER          , INTENT(IN)    :: k_tend
    REAL             , INTENT(IN)    :: rho_air(m1,m2,m3)

    INTEGER i,j

    !- also this routine assumes that the source term is expressed in terms of  density (kg/m3)
    !- and, then, the conversion to mixing ratio (division by air density) is needed.
    !- k_src express the level where the src data is stored in the sc_src array
    ! if(nint(alfa) == 0) then

    !   do j=ja,jz ; do i=ia,iz
    !  
    !
    !   sc_t(k_tend,i,j) = sc_t(k_tend,i,j) +  (                        &
    !                      sc_src1(k_src,i,j)* ( 1.0D0- rt(i,j) )*alfa  +  &
    !                      sc_src2(k_src,i,j)* rt(i,j)               )/ &
    !                      rho_air(k_tend,i,j) !- the level of air density
    !	 	                         !- must be the same of tendency
    !                     
    !   enddo;enddo

    !original
    IF(alfa > 0.) THEN 


       sc_t(k_tend,ia:iz,ja:jz) = sc_t   (k_tend ,ia:iz,ja:jz) + (  &
            sc_src1(k_src  ,ia:iz,ja:jz) * ( 1.0D0- rt(ia:iz,ja:jz) )*alfa +     &
            sc_src2(k_src  ,ia:iz,ja:jz) *          rt(ia:iz,ja:jz)        )  /  &
            rho_air(k_tend ,ia:iz,ja:jz) !- the level of air density
       !- must be the same of tendency
    ELSE

       sc_t(k_tend,ia:iz,ja:jz) = sc_t   (k_tend ,ia:iz,ja:jz) + (  &
!!!                           sc_src1(k_src  ,ia:iz,ja:jz) * ( 1.0D0- rt(ia:iz,ja:jz) )*alfa  +  &
            sc_src2(k_src  ,ia:iz,ja:jz) *          REAL(rt(ia:iz,ja:jz))         )  /  &
            rho_air(k_tend ,ia:iz,ja:jz) !- the level of air density
       !- must be the same of tendency

    ENDIF

  END SUBROUTINE source_to_tend_cycle

  !------------------------------------------------------------------
  ! NOT USED
  !
  !  subroutine invert_air_dens(m1,rho_air)
  !   implicit none
  !   integer m1
  !   real rho_air(m1)
  !   rho_air(:)=1./rho_air(:)
  ! 
  !  end subroutine invert_air_dens
  !------------------------------------------------------------------


  SUBROUTINE alloc_emiss_cycle(mmxp,mmyp,ngrids,nsrc,maxgrds)

    ! node_mod
    INTEGER , INTENT(IN) :: mmxp(maxgrds)
    INTEGER , INTENT(IN) :: mmyp(maxgrds)
    
    ! mem_grid
    INTEGER , INTENT(IN) :: ngrids
    
    ! mem_chem1
    INTEGER , INTENT(IN) :: nsrc
    
    ! grid_dims
    INTEGER, INTENT(IN) :: maxgrds

    INTEGER ng,isrc

    ALLOCATE (emiss_cycle(nsrc,ngrids))
    DO ng=1,ngrids
       DO isrc=1,nsrc
          ALLOCATE(emiss_cycle(isrc,ng)%dcnorma_inv( mmxp(ng),mmyp(ng) ) )    
          emiss_cycle(isrc,ng)%dcnorma_inv = 0.
          ALLOCATE(emiss_cycle(isrc,ng)%emission_rate( mmxp(ng),mmyp(ng) ) )    
          emiss_cycle(isrc,ng)%emission_rate = 0.

       ENDDO
    ENDDO

    emiss_cycle_alloc=.TRUE.

  END SUBROUTINE alloc_emiss_cycle

  !------------------------------------------------------------------
  SUBROUTINE get_diurnal_cycle_normalized(ngrid,m2,m3,ia,iz,ja,jz,dtlt,glat, &
                                          glon,imonth1,idate1,iyear1,itime1, &
                                          antro,bioge,pi180)
      
    ! original
    INTEGER , INTENT(IN) :: ngrid
    INTEGER , INTENT(IN) :: m2
    INTEGER , INTENT(IN) :: m3
    INTEGER , INTENT(IN) :: ia
    INTEGER , INTENT(IN) :: iz
    INTEGER , INTENT(IN) :: ja
    INTEGER , INTENT(IN) :: jz
    REAL    , INTENT(IN) :: dtlt
    REAL    , INTENT(IN) :: glat(m2,m3)
    REAL    , INTENT(IN) :: glon(m2,m3)
    
    ! mem_grid
    INTEGER , INTENT(IN) :: imonth1
    INTEGER , INTENT(IN) :: idate1
    INTEGER , INTENT(IN) :: iyear1
    INTEGER , INTENT(IN) :: itime1
      
    ! mem_chem1
    INTEGER , INTENT(IN) :: antro
    INTEGER , INTENT(IN) :: bioge

    ! rconstants
    REAL    , INTENT(IN) :: pi180

    INTEGER :: jday,i,j,julday

    REAL :: solfac,tdec,sdec,cdec,declin,d0,d02,dayhr,radlat,cslcsd,snlsnd, &
            gglon,dayhrr,hrangl

    !- local var
    INTEGER it
    REAL time_x,dt_x,cosz
    REAL dcnorma(m2,m3)
    REAL xxx

    !print*,'-----------------------------------------------------' 
    !print*,'getting diurnal/space variable emission rates - time=',time/3600. ;call flush(6)
    !print*,'-----------------------------------------------------'


    dt_x=0.
    time_x=0.
    dcnorma(:,:)=0.


    DO it=1,NINT(86400./dtlt) 

       jday = julday(imonth1,idate1,iyear1) ! don't change the position of this line

       jday = jday + NINT(time_x/86400.)
       !      sdec - sine of declination, cdec - cosine of declination
       declin = -23.5 * COS(6.283 / 365. * (jday + 9)) * pi180
       sdec = SIN(declin)
       cdec = COS(declin)

       ! Find the factor, solfac, to multiply the solar constant to correct
       ! for Earth's varying distance to the sun.

       !d0 = 6.2831853 * float(jday-1) / 365.
       !d02 = d0 * 2.
       !solfac = 1.000110 + 0.034221 * cos (d0) + 0.001280 * sin(d0)  &
       !     + 0.000719 * cos(d02) + 0.000077 * sin(d02)

       ! Find the hour angle, THEN get cosine of zenith angle.

       dayhr = time_x / 3600. + float(itime1/100) + float(MOD(itime1,100)) / 60.

       DO j = ja,jz
          DO i = ia,iz
             radlat = glat(i,j) * pi180
             !IF (lonrad .eq. 0) radlat = centlat(1) * pi180
             !IF (radlat .eq. declin) radlat = radlat + 1.e-5
             cslcsd = COS(radlat) * cdec
             snlsnd = SIN(radlat) * sdec
             !gglon = glon(i,j)
             !IF (lonrad .eq. 0) gglon = centlon(1)
             dayhrr = MOD(dayhr+glon(i,j)/15.+24.,24.)
             hrangl = 15. * (dayhrr - 12.) * pi180
             cosz = snlsnd + cslcsd * COS(hrangl)
             !cosz = min(cosz+1.0E-10, 1.0) 

             dcnorma(i,j)=dcnorma(i,j)+MAX(0.,cosz)

          END DO
       END DO
       time_x=time_x+dtlt
    END DO

    !- invert dcnorma to save computation time
    dcnorma(:,:)=1./(dcnorma(:,:)*dtlt)
    !- transfer the emission cycle to bioge and antro arrays only, for now  
    emiss_cycle(bioge,ngrid)%dcnorma_inv( :,:) = dcnorma(:,:)
    emiss_cycle(antro,ngrid)%dcnorma_inv( :,:) = dcnorma(:,:)

    !print*,'dcnorma_INV=',dcnorma(3,4); call flush(6)

  END SUBROUTINE get_diurnal_cycle_normalized


  !------------------------------------------------------------------
  ! (DMK) NOT USED
  !
  ! SUBROUTINE get_cosz(m2,m3,ia,iz,ja,jz,glat,glon,cosz)
  !
  !  use mem_grid   , only: imonth1,idate1,iyear1,time,itime1,centlat, &
  !       centlon
  !  use mem_radiate, only: lonrad
  !  use rconstants , only: pi180
  !
  !  implicit none
  !
  !  integer :: m2,m3,ia,iz,ja,jz,jday,i,j,julday
  !
  !  real :: solfac,tdec,sdec,cdec,declin,d0,d02,dayhr,radlat,cslcsd,snlsnd  &
  !       ,gglon,dayhrr,hrangl,dtlt
  !  real, dimension(m2,m3) :: glat,glon,cosz
  !
  !  jday = julday(imonth1,idate1,iyear1)
  !  
  !  jday = jday + nint(time/86400.)
  !  !      sdec - sine of declination, cdec - cosine of declination
  !  declin = -23.5 * cos(6.283 / 365. * (jday + 9)) * pi180
  !  sdec = sin(declin)
  !  cdec = cos(declin)
  !  
  !  ! Find the factor, solfac, to multiply the solar constant to correct
  !  ! for Earth's varying distance to the sun.
  !  
  !  !d0 = 6.2831853 * float(jday-1) / 365.
  !  !d02 = d0 * 2.
  !  !solfac = 1.000110 + 0.034221 * cos (d0) + 0.001280 * sin(d0)  &
  !  !     + 0.000719 * cos(d02) + 0.000077 * sin(d02)
  !  
  !  ! Find the hour angle, THEN get cosine of zenith angle.
  !  
  !  dayhr = time / 3600. + float(itime1/100) + float(mod(itime1,100)) / 60.
  !  
  !  DO j = ja,jz
  !     DO i = ia,iz
  !        radlat = glat(i,j) * pi180
  !        !IF (lonrad .eq. 0) radlat = centlat(1) * pi180
  !        !IF (radlat .eq. declin) radlat = radlat + 1.e-5
  !        cslcsd = cos(radlat) * cdec
  !        snlsnd = sin(radlat) * sdec
  !        !gglon = glon(i,j)
  !        !IF (lonrad .eq. 0) gglon = centlon(1)
  !        dayhrr = mod(dayhr+glon(i,j)/15.+24.,24.)
  !        hrangl = 15. * (dayhrr - 12.) * pi180
  !        cosz(i,j) = snlsnd + cslcsd * cos(hrangl)
  !        
  !     END DO
  !  END DO
  !  
  ! END SUBROUTINE get_cosz

  !------------------------------------------------------------------
  ! (DMK) NOT USED
  ! SUBROUTINE Greg2Jul(h, d, m, y, jd)
  !  
  !   INTEGER , INTENT(IN)  :: d
  !   INTEGER , INTENT(IN)  :: m
  !   INTEGER , INTENT(IN)  :: y
  !   INTEGER , INTENT(IN)  :: h
  !   REAL    , INTENT(OUT) :: jd
  !
  !   jd = ( 1461 * ( y + 4800 + ( m - 14 ) / 12 ) ) / 4 +        &
  !        ( 367 * ( m - 2 - 12 * ( ( m - 14 ) / 12 ) ) ) / 12 -  &
  !        ( 3 * ( ( y + 4900 + ( m - 14 ) / 12 ) / 100 ) ) / 4 + &
  !        d - 32075
  !  
  !   jd = jd + (h/24.)
  !  
  ! END SUBROUTINE Greg2Jul
  !------------------------------------------------------------------

END MODULE chem_sources




























!----------------------------------------------------------------------------
! (DMK) NOT USED
!
! !- to be used for interpolation from some data (lat/lon gridded) to model 
! SUBROUTINE interp_init_concentration(ng,n1,n2,n3,glat,glon,dn0,rtgt,dxt,dyt)
!  IMPLICIT NONE
!
!  INTEGER,INTENT(IN) :: ng,n1,n2,n3
!
!  REAL, DIMENSION(n2,n3)    :: dxt,dyt,rtgt,glat,glon
!  REAL, DIMENSION(n1,n2,n3) :: s3p,dn0
!  !variaveis locais
!  REAL ::  latnf,lonnf,dlatr,dlonr  !,ilats,lonn,lats,ilons,latn,lons
!  REAL ::  undef  
!  INTEGER :: i,j,k,qi1,qi2,qj1,qj2,ncount,ii,jj,jc,ic,i1,j1,i2,j2,nzvar !,kk
!  REAL, ALLOCATABLE :: DATAf(:,:)	 ! model concentration
!  REAL, ALLOCATABLE :: tmp(:,:,:)	 !dado dummy temporario
!  !- informacoes da grade dos dados que serao interpolados (input data)
!  REAL, ALLOCATABLE :: api_us(:,:,:),prlat(:,:),prlon(:,:),usdum(:) !api_us= dado inicial a ser interpolado
!  INTEGER,PARAMETER ::nlon=546,nlat=596,n4us=1  !n4us eh usado no caso do dado nao ser somente de um nivel vertical
!  REAL,PARAMETER :: ilatn= 0.090497! delta lat
!  REAL,PARAMETER :: ilonn= 0.091743! delta lon
!  REAL,PARAMETER :: latni=-38.! lat mais ao sul
!  REAL,PARAMETER :: lonni=-85.! lon mais a oeste
!  !real,parameter :: fcu =1.e+6  !=> mg [gas/part] /kg [ar]
!
!  latnf=latni + (nlat-1)*ilatn
!  lonnf=lonni + (nlon-1)*ilonn
!
!  !valor indefinido para o modelo
!  undef = 0.
!
!
!  !grade do dado a ser interpolado
!  ALLOCATE(prlat(nlon,nlat),prlon(nlon,nlat))
!  CALL api_prlatlon(nlon,nlat,prlat,prlon,ilatn,ilonn,latni,lonni)
!  ALLOCATE(api_us(nlon,nlat,n4us),usdum(n4us))
!
!  !- need to adapt the lines below in order to read the data:
!  PRINT*,'--------------------------------------'
!  PRINT*,'Opening initial data='
!  nzvar=22 !number of vert levels of data to be ingested
!  ALLOCATE(tmp(nlon,nlat,nzvar))
!
!  OPEN(2,status='OLD',form='unformatted',access='direct', &
!       recl=4*nlat*nlon*nzvar,file='9_19_2002_262.bin')
!  READ(2) tmp
!  CLOSE(2)
!  !	   call swap32(tmp,nlat*nlon*nzvar)
!  DO i=1,nlon
!     DO j=1,nlat
!	api_us(i,j,1) = tmp (i,j,18)
!	!	   if(tmp(i,j,18).gt. 0.)print*,i,j,tmp (i,j,18)
!     ENDDO
!  ENDDO
!  PRINT*,'--------------------------------------'
!
!  ! dado final        na grade do modelo
!  ALLOCATE(DATAf(n2,n3))
!
!
!  ! loop no dominio do modelo
!  DO j = 1,n3
!     DO i = 1,n2
!
!	! evite pontos fora do dominio da grade de umidade de solo
!	IF(glat(i,j) .LT. latni .OR. glat(i,j) .GT. latnf .OR. &
!	     glon(i,j) .LT. lonni .OR. glon(i,j) .GT. lonnf) THEN
!
!	   DATAf(i,j) = undef
!	   GOTO 1111
!	ENDIF
!
!	CALL interpolacao (glon(i,j),glat(i,j),nlon,nlat,prlat,prlon,i1,i2,ic,j1,j2,jc)
!	IF(ic.GE.0 .AND. jc .GE. 0) THEN
!	   !print*,ic,jc,i,j,ifm
!	   dlonr=0.5*(glon(n2,j)-glon(1,j))/float(n2-1)
!	   dlatr=0.5*(glat(i,n3)-glat(i,1))/float(n3-1)
!	   qi1=INT(dlonr/ilonn+0.5)
!	   qi2=INT(dlonr/ilonn+0.5)
!	   qj1=INT(dlatr/ilatn+0.5)
!	   qj2=INT(dlatr/ilatn+0.5)
!	   !print*,jc-qj1,jc+qj2,ic-qi1,ic+qi2
!	   DO k=1,n4us
!	      ncount = 0
!	      usdum(k)=0.
!
!	      ! print*,'======================================== i j=',i,j
!	      ! print*,jc-qj1,jc+qj2,ic-qi1,ic+qi2
!	      DO jj =MAX(1,jc-qj1),MIN(nlat,jc+qj2)
!		 DO ii = MAX(1,ic-qi1),MIN(nlon,ic+qi2)
!
!		    IF (api_us(ii,jj,k).GE.1.e-5) THEN
!		       ncount = ncount + 1
!		       usdum(k) = usdum(k) + api_us(ii,jj,k)
!
!		       ! print*,'ii jj=',ii,jj
!		       ! print*,'ncount datai soma=',ncount,api_us(ii,jj,k),usdum(k)
!
!		    ENDIF
!		 ENDDO
!	      ENDDO
!	      IF(ncount .GT. 0 ) THEN
!		 DATAf(i,j) = usdum(k) / (float(ncount))
!	      ELSE
!		 DATAf(i,j) = undef
!	      ENDIF
!
!	      ! print*,'FINAL=',dataf(i,j)
!	      ! print*,'========================================'
!
!	      !if(usdum(k) .lt. 1.e-5) then
!	      ! print*,glon(i,j),glat(i,j),usdum(k)
!	      ! stop
!	      !endif
!	   ENDDO
!	   !
!	ENDIF
! 1111	CONTINUE
!
!     ENDDO
!  ENDDO
!
!  DEALLOCATE(api_us,usdum,prlat,prlon)
!
!  OPEN(2,status='UNKNOWN',form='unformatted',access='direct', &
!       recl=4*n2*n3,file='new.bin')
!  WRITE(2) DATAf
!  CLOSE(2)
!
!  
!  CALL azero(n1*n2*n3,s3p)   !zera s3p
!  DO j=1,n3
!     DO i=1,n2
!	DO k=1,7
!	   s3p(k,i,j)  = DATAf(i,j)/dn0(k,i,j)   
!	ENDDO
!     ENDDO
!  ENDDO
!  !
!  DEALLOCATE(tmp,DATAf)
! END SUBROUTINE interp_init_concentration
!-----------------------------------------------------------------------
